Two-way interaction between arboviruses and the microbiota of mosquitoes
1 Eukaryotic virome analysis
#### load phyloseq project ####
GP.M1 <- readRDS("~/Desktop/mosq_Guad_infection/RDS/GIM_virome_raw.RDS")
GP.M1## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 152 taxa and 154 samples ]
## sample_data() Sample Data: [ 154 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 152 taxa by 7 taxonomic ranks ]1.1 Eukaryotic virome of Aedes aegypti
#### color ####
library(viridis)
col = plasma(9,alpha = 1, begin=0.9, end = 0, direction = 1)
scales::show_col(col)
1.1.1 Subset eukaryotic virome of Aedes aegypti
GP.M1_AA <- subset_samples(GP.M1, mosq_species == "Ae_aegypti")
GP.M1_AA <- subset_samples(GP.M1_AA, group.plaque.qPCR.1000. != "")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Chuvirus Mos8Chu0")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Kaiowa virus")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Cumbaru virus")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Guato virus")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Trichoplusia ni TED virus")
GP.M1_AA <- subset_taxa(GP.M1_AA, species != "Guadeloupe Culex tymo-like virus")
GP.M1_AA## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 126 taxa and 92 samples ]
## sample_data() Sample Data: [ 92 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 126 taxa by 7 taxonomic ranks ]sample_data(GP.M1_AA)$Treatment <- factor(sample_data(GP.M1_AA)$Treatment, levels = c("AA-sucrose", "AA-blood", "AA-ZIKV"))
sample_data(GP.M1_AA)$Info1 <- factor(sample_data(GP.M1_AA)$Info1, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV", "AA-21dpi-sucrose","AA-21dpi-blood","AA-21dpi-ZIKV"))
sample_data(GP.M1_AA)$group.plaque.qPCR.0.body.0. <- factor(sample_data(GP.M1_AA)$group.plaque.qPCR.0.body.0, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV(-/-/-)", "AA-7dpi-ZIKV(-/-/+)", "AA-7dpi-ZIKV(-/+/+)", "AA-21dpi-sucrose", "AA-21dpi-blood", "AA-21dpi-ZIKV(-/-/-)", "AA-21dpi-ZIKV(-/-/+)", "AA-21dpi-ZIKV(-/+/-)", "AA-21dpi-ZIKV(-/+/+)", "AA-21dpi-ZIKV(+/+/+)"))
sample_data(GP.M1_AA)$plaque_assay <- factor(sample_data(GP.M1_AA)$plaque_assay,
levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV(-)",
"AA-21dpi-sucrose", "AA-21dpi-blood", "AA-21dpi-ZIKV(-)", "AA-21dpi-ZIKV(+)"))
sample_data(GP.M1_AA)$dpi <- factor(sample_data(GP.M1_AA)$dpi, levels = c("7dpi", "21dpi"))1.1.2 Merge on species level
GP.M1_species <- GP.M1_AA %>% tax_glom(taxrank = "species")
GP.M1_species## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 55 taxa and 92 samples ]
## sample_data() Sample Data: [ 92 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 55 taxa by 7 taxonomic ranks ]GP.M1_samselect = prune_samples(sample_sums(GP.M1_species) >0, GP.M1_species)
# species
GP.M1_species_select_0 = prune_taxa(taxa_sums(GP.M1_samselect) > 0, GP.M1_samselect)
GP.M1_species_select_0## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 40 taxa and 92 samples ]
## sample_data() Sample Data: [ 92 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 40 taxa by 7 taxonomic ranks ]GP.M1_species_select = prune_taxa(taxa_sums(GP.M1_samselect) > 500, GP.M1_samselect)
GP.M1_species_select## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 19 taxa and 92 samples ]
## sample_data() Sample Data: [ 92 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 19 taxa by 7 taxonomic ranks ]1.1.3 Heatmap on viral species level
#### get table of reads for each viral species ####
OTU_species <- merge(otu_table(GP.M1_species_select), tax_table(GP.M1_species_select), by=0)
rownames(OTU_species) <- OTU_species$species
#### order viral species by familiy ####
OTU_species$family <- as.character(OTU_species$family)
OTU_species$family## [1] "Virgaviridae" "Rhabdoviridae" "Flaviviridae" "Orthomyxoviridae"
## [5] "Flaviviridae" "Totiviridae" NA "Xinmoviridae"
## [9] "Totiviridae" NA "Totiviridae" "Flaviviridae"
## [13] NA "Reoviridae" "Flaviviridae" NA
## [17] "Phenuiviridae" "Flaviviridae" "Xinmoviridae"OTU_species$family[is.na(OTU_species$family)] <- "unclassified"
OTU_species_m<- melt(OTU_species)## Using Row.names, superkingdom, phylum, class, order, family, genus, species as id variablesOTU_species_m_agg <- aggregate(value~family+species, FUN=sum, data=OTU_species_m)
OTU_species_m_agg1 <- aggregate(value~family, FUN=sum, data=OTU_species_m)
family_ord <- OTU_species_m_agg1[order(-OTU_species_m_agg1$value),]$family
family_ord <- c(family_ord[-3],"unclassified")
family_ord## [1] "Phenuiviridae" "Totiviridae" "Flaviviridae" "Orthomyxoviridae"
## [5] "Xinmoviridae" "Virgaviridae" "Reoviridae" "Rhabdoviridae"
## [9] "unclassified"OTU_species_m_agg$index <- 1
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Totiviridae")] <- 2
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Flaviviridae")] <- 3
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Orthomyxoviridae")] <- 4
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Xinmoviridae")] <- 5
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Metaviridae")] <- 6
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Virgaviridae")] <- 7
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Rhabdoviridae")] <- 8
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Reoviridae")] <- 9
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "unclassified")] <- 10
OTU_species_m_agg <- OTU_species_m_agg[order(OTU_species_m_agg$index, -OTU_species_m_agg$value),]
OTU_species_m_agg## family species value index
## 14 Phenuiviridae Phasi Charoen-like phasivirus 8253667 1
## 2 Totiviridae Aedes aegypti totivirus 7849526 2
## 1 Totiviridae Aedes aegypti toti-like virus 43484 2
## 5 Totiviridae Australian Anopheles totivirus 1722 2
## 19 Flaviviridae Zika virus 409642 3
## 8 Flaviviridae Cell fusing agent virus 55036 3
## 16 Flaviviridae unclassified Flaviviridae 8111 3
## 13 Flaviviridae Menghai flavivirus 6987 3
## 18 Flaviviridae Xishuangbanna aedes flavivirus 915 3
## 9 Orthomyxoviridae Guadeloupe mosquito quaranja-like virus 1 57556 4
## 4 Xinmoviridae Aedes anphevirus 21985 5
## 15 Xinmoviridae unclassified Anphevirus 18618 5
## 7 Virgaviridae Bombus-associated virus Vir1 1644 7
## 17 Rhabdoviridae unclassified Rhabdoviridae 741 8
## 12 Reoviridae Liao ning virus 755 9
## 10 unclassified Guadeloupe mosquito virus 583366 10
## 6 unclassified Beihai barnacle virus 12 231777 10
## 11 unclassified Humaita-Tubiacanga virus 7435 10
## 3 unclassified Aedes alboannulatus orthomyxo-like virus 868 10species_ord <- as.character(OTU_species_m_agg$species)
species_ord## [1] "Phasi Charoen-like phasivirus"
## [2] "Aedes aegypti totivirus"
## [3] "Aedes aegypti toti-like virus"
## [4] "Australian Anopheles totivirus"
## [5] "Zika virus"
## [6] "Cell fusing agent virus"
## [7] "unclassified Flaviviridae"
## [8] "Menghai flavivirus"
## [9] "Xishuangbanna aedes flavivirus"
## [10] "Guadeloupe mosquito quaranja-like virus 1"
## [11] "Aedes anphevirus"
## [12] "unclassified Anphevirus"
## [13] "Bombus-associated virus Vir1"
## [14] "unclassified Rhabdoviridae"
## [15] "Liao ning virus"
## [16] "Guadeloupe mosquito virus"
## [17] "Beihai barnacle virus 12"
## [18] "Humaita-Tubiacanga virus"
## [19] "Aedes alboannulatus orthomyxo-like virus"which(species_ord == "Zika virus")## [1] 5species_ord_1 <- c("Zika virus", species_ord[-5])
#### transform the table - reads for each viral species #####
rm <- c("Row.names", "superkingdom", "phylum", "class", "order", "family", "genus", "species")
OTU_species <- OTU_species[,!colnames(OTU_species) %in% rm]
OTU_species_trans <- log10(OTU_species)
is.na(OTU_species_trans) <- sapply(OTU_species_trans, is.infinite)
OTU_species_trans[is.na(OTU_species_trans)]<-0
df.original <- OTU_species
df <- OTU_species_trans
#### use meta data for top annotation ####
meta = read.csv("~/Desktop/mosq_Guad_infection/GIM_metadata_new.csv", header=TRUE, row.names=1, sep=",")
meta_s = meta[!meta$mosq_species == "",]
meta_AA <- meta[colnames(df),]
#### order sample by ZIKV reads within each group & column annotation ####
meta_AAz <- merge(meta_AA, as.data.frame(t(df)),by=0)
head(meta_AAz[1:5, 1:5])## Row.names mosq body mosq_or head
## 1 GIM100 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 2 GIM22 AaZ-7dpi-2 GIM-22 AMT-7dpi-2 GIM-2
## 3 GIM23 AaZ-7dpi-3 GIM-23 AMT-7dpi-3 GIM-3
## 4 GIM24 AaZ-7dpi-4 GIM-24 AMT-7dpi-4 GIM-4
## 5 GIM25 AaZ-7dpi-5 GIM-25 AMT-7dpi-5 GIM-5meta_AAz$Treatment <- as.character(meta_AAz$Treatment)
meta_AAz$Treatment <- factor(meta_AAz$Treatment, levels = c("AA-sucrose", "AA-blood", "AA-ZIKV"))
meta_AAz$Info1 <- as.character(meta_AAz$Info1)
meta_AAz$Info1 <- factor(meta_AAz$Info1, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV", "AA-21dpi-sucrose","AA-21dpi-blood","AA-21dpi-ZIKV"))
meta_AAz$qPCR_head._ZIKV_WNV <- log10(meta_AAz$qPCR_head._ZIKV_WNV)
meta_AAz$qPCR_body._ZIKV_WNV <- log10(meta_AAz$qPCR_body._ZIKV_WNV)
which(colnames(meta_AAz)=="qPCR_head._ZIKV_WNV")## [1] 16which(colnames(meta_AAz)=="qPCR_body._ZIKV_WNV")## [1] 17is.na(meta_AAz[,16:17])<-sapply(meta_AAz[,16:17], is.infinite)
meta_AAz[,16:17][is.na(meta_AAz[,16:17])]<-0
meta_AAz$group_n1 <- meta_AAz$plaque_assay
meta_AAz$plaque_assay <- factor(meta_AAz$plaque_assay, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV(-)", "AA-21dpi-sucrose", "AA-21dpi-blood", "AA-21dpi-ZIKV(-)","AA-21dpi-ZIKV(+)"), ordered = TRUE)
meta_AAz_order1 <- meta_AAz[meta_AAz$dpi== "7dpi",][with(meta_AAz[meta_AAz$dpi== "7dpi",], order(plaque_assay, -qPCR_head._ZIKV_WNV, -qPCR_body._ZIKV_WNV, -`Zika virus`, -`Guadeloupe mosquito virus`,-`Phasi Charoen-like phasivirus`, -`Aedes aegypti toti-like virus`,-`Aedes anphevirus`,-`Guadeloupe mosquito quaranja-like virus 1`)), ]
meta_AAz_order2 <- meta_AAz[meta_AAz$dpi== "21dpi",][with(meta_AAz[meta_AAz$dpi== "21dpi",], order(plaque_assay, -qPCR_head._ZIKV_WNV, -qPCR_body._ZIKV_WNV, -`Zika virus`, -`Guadeloupe mosquito virus`,-`Phasi Charoen-like phasivirus`,-`Aedes aegypti toti-like virus`,-`Aedes anphevirus`,-`Guadeloupe mosquito quaranja-like virus 1`)), ]
meta_AAz_order <- rbind(meta_AAz_order1,meta_AAz_order2)
df <- df[,meta_AAz_order$Row.names]
col_ha = HeatmapAnnotation(qPCR = anno_lines(meta_AAz_order[,16:17], gp = gpar(col = c("black","darkgrey")), height = unit(1.5, "cm"),ylim = c(-0.5,8), axis_param = list(at = c(0,2,4,6,8)),add_points = TRUE, pt_gp = gpar(col = c("black","darkgrey")), pch = c(16, 17)))
#### column split ####
df_1 <- df
df_1$ZIKV <- "others"
df_1$ZIKV[which(rownames(df_1)=="Zika virus")] <- "aZIKV"#### draw figure #####
heatmap_col_tp2 <- colorRampPalette(brewer.pal(9, "YlGnBu"), space = "rgb")(100)
Heatmap(as.matrix(df),
name = "log10 reads number", #title of legend
column_title = "Samples", row_title = "viral species",
row_names_gp = gpar(fontsize = 10),
column_names_gp = gpar(fontsize = 6),
column_title_gp = gpar(fontsize = 12),
cluster_rows = F,
cluster_columns= F,
row_order = species_ord_1,
column_labels = meta_AAz_order$mosq,
top_annotation = col_ha,
col = heatmap_col_tp2,
column_split = meta_AAz_order$group_n1,
row_split = df_1$ZIKV,
border = TRUE,
rect_gp = gpar(col = "white", lwd = 0.5))
1.1.4 Alpha diversity on species level
plot_richness <- function (physeq, x = "samples", color = NULL, shape = NULL,
title = NULL, scales = "free_y", nrow = 1, shsi = NULL, measures = NULL,
sortby = NULL)
{
erDF = estimate_richness(physeq, split = TRUE, measures = measures)
measures = colnames(erDF)
ses = colnames(erDF)[grep("^se\\.", colnames(erDF))]
measures = measures[!measures %in% ses]
if (!is.null(sample_data(physeq, errorIfNULL = FALSE))) {
DF <- data.frame(erDF, sample_data(physeq))
}
else {
DF <- data.frame(erDF)
}
if (!"samples" %in% colnames(DF)) {
DF$samples <- sample_names(physeq)
}
if (!is.null(x)) {
if (x %in% c("sample", "samples", "sample_names", "sample.names")) {
x <- "samples"
}
}
else {
x <- "samples"
}
mdf = reshape2::melt(DF, measure.vars = measures)
mdf$se <- NA_integer_
if (length(ses) > 0) {
selabs = ses
names(selabs) <- substr(selabs, 4, 100)
substr(names(selabs), 1, 1) <- toupper(substr(names(selabs),
1, 1))
mdf$wse <- sapply(as.character(mdf$variable), function(i,
selabs) {
selabs[i]
}, selabs)
for (i in 1:nrow(mdf)) {
if (!is.na(mdf[i, "wse"])) {
mdf[i, "se"] <- mdf[i, (mdf[i, "wse"])]
}
}
mdf <- mdf[, -which(colnames(mdf) %in% c(selabs, "wse"))]
}
if (!is.null(measures)) {
if (any(measures %in% as.character(mdf$variable))) {
mdf <- mdf[as.character(mdf$variable) %in% measures,
]
}
else {
warning("Argument to `measures` not supported. All alpha-diversity measures (should be) included in plot.")
}
}
if (!is.null(shsi)) {
warning("shsi no longer supported option in plot_richness. Please use `measures` instead")
}
if (!is.null(sortby)) {
if (!all(sortby %in% levels(mdf$variable))) {
warning("`sortby` argument not among `measures`. Ignored.")
}
if (!is.discrete(mdf[, x])) {
warning("`sortby` argument provided, but `x` not a discrete variable. `sortby` is ignored.")
}
if (all(sortby %in% levels(mdf$variable)) & is.discrete(mdf[,
x])) {
wh.sortby = which(mdf$variable %in% sortby)
mdf[, x] <- factor(mdf[, x], levels = names(sort(tapply(X = mdf[wh.sortby,
"value"], INDEX = mdf[wh.sortby, x], mean, na.rm = TRUE,
simplify = TRUE))))
}
}
richness_map = aes_string(x = x, y = "value", colour = color,
shape = shape)
p = ggplot(mdf, richness_map) + geom_boxplot(na.rm = TRUE,outlier.shape = NA)+
# geom_point(size=1)
geom_jitter(position=position_jitter(0.25), size = 3.5, alpha=0.7)
# if (any(!is.na(mdf[, "se"]))) {
# p = p + geom_errorbar(aes(ymax = value + se, ymin = value -
# se), width = 0.1)
# }
# p = p + theme(axis.text.x = element_text(angle = -90, vjust = 0.5,
# hjust = 0))
p = p + ylab("Alpha Diversity Measure")
p = p + facet_wrap(~variable, nrow = nrow, scales = scales)
if (!is.null(title)) {
p <- p + ggtitle(title)
}
return(p)
}
#### on diet ####
plot_richness(GP.M1_species_select_0, x="Info1",measures=c("Shannon"),col="Treatment", shape = "dpi") +
scale_color_manual(values = rev(c(col[9],col[5],col[2]))) +
scale_shape_manual(values = c(16, 17))+
ylim(-0.001, 2.5)+
theme_bw()+
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 12, colour = "black", angle = 90, hjust = 1),
strip.text.x = element_text(size = 12, colour = "black"),
plot.title = element_text(hjust = 0.5, size = 12)) +
ggtitle("Alpha diversity of eukaryotic virus in Ae. aegypti")+
ylab("Alpha Diversity Value")
erich <- estimate_richness(GP.M1_species_select_0, measures = c("Shannon"))
pairwise.wilcox.test(erich$Shannon,sample_data(GP.M1_species_select_0)$Info1, p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: erich$Shannon and sample_data(GP.M1_species_select_0)$Info1
##
## AA-7dpi-sucrose AA-7dpi-blood AA-7dpi-ZIKV AA-21dpi-sucrose
## AA-7dpi-blood 0.57359 - - -
## AA-7dpi-ZIKV 4.1e-05 0.03450 - -
## AA-21dpi-sucrose 0.02381 1.00000 0.00437 -
## AA-21dpi-blood 0.03450 0.68452 0.24051 0.73980
## AA-21dpi-ZIKV 0.00092 0.57359 4.1e-05 0.41492
## AA-21dpi-blood
## AA-7dpi-blood -
## AA-7dpi-ZIKV -
## AA-21dpi-sucrose -
## AA-21dpi-blood -
## AA-21dpi-ZIKV 0.70485
##
## P value adjustment method: BH#### on diet and infection status ####
sample_data(GP.M1_species_select_0)$hd <- as.character(sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000.)
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-7dpi-ZIKV(-/+/+)"] <- "7dpe-Head(+)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-7dpi-ZIKV(-/-/+)"] <- "7dpe-Head(-)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-7dpi-ZIKV(-/-/-)"] <- "7dpe-Head(-)/Body(-)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-21dpi-ZIKV(+/+/+)"] <- "21dpe-Head(+)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-21dpi-ZIKV(-/+/+)"] <- "21dpe-Head(+)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-21dpi-ZIKV(-/-/+)"] <- "21dpe-Head(-)/Body(+)"
sample_data(GP.M1_species_select_0)$hd[sample_data(GP.M1_species_select_0)$group.plaque.qPCR.1000.body.1000. == "AA-21dpi-ZIKV(-/-/-)"] <- "21dpe-Head(-)/Body(-)"
sample_data(GP.M1_species_select_0)$hd <- factor(sample_data(GP.M1_species_select_0)$hd,
levels = c("AA-7dpi-sucrose","AA-7dpi-blood","7dpe-Head(-)/Body(-)","7dpe-Head(-)/Body(+)", "7dpe-Head(+)/Body(+)",
"AA-21dpi-sucrose","AA-21dpi-blood","21dpe-Head(-)/Body(-)","21dpe-Head(-)/Body(+)", "21dpe-Head(+)/Body(+)"))
plot_richness(GP.M1_species_select_0, x="hd", measures=c("Shannon"), col="hd", shape = "dpi") +
scale_color_manual(values = rev(c(col[9],col[7],col[5],col[3],col[1],col[9],col[7],col[5],col[3],col[1]))) +
scale_shape_manual(values = c(16, 17))+
ylim(-0.0001,2.5)+
theme_bw()+
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 12, colour = "black", angle = 90, hjust = 1, vjust=0.5),
strip.text.x = element_text(size = 12, colour = "black"),
plot.title = element_text(hjust = 0.5, size = 12)) +
ggtitle("Alpha diversity of eukaryotic virus in Ae. aegypti")+
ylab("Alpha Diversity Value")
erich <- estimate_richness(GP.M1_species_select_0, measures = c("Shannon"))
pairwise.wilcox.test(erich$Shannon,sample_data(GP.M1_species_select_0)$hd, p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: erich$Shannon and sample_data(GP.M1_species_select_0)$hd
##
## AA-7dpi-sucrose AA-7dpi-blood 7dpe-Head(-)/Body(-)
## AA-7dpi-blood 0.6760 - -
## 7dpe-Head(-)/Body(-) 0.0013 0.0887 -
## 7dpe-Head(-)/Body(+) 0.0112 0.1111 0.9481
## 7dpe-Head(+)/Body(+) 0.0549 0.2500 0.7858
## AA-21dpi-sucrose 0.0397 1.0000 0.0199
## AA-21dpi-blood 0.0549 0.7766 0.2983
## 21dpe-Head(-)/Body(-) 0.0199 0.6814 0.0149
## 21dpe-Head(-)/Body(+) 0.5669 0.7895 1.0000
## 21dpe-Head(+)/Body(+) 0.0013 0.6814 0.0054
## 7dpe-Head(-)/Body(+) 7dpe-Head(+)/Body(+)
## AA-7dpi-blood - -
## 7dpe-Head(-)/Body(-) - -
## 7dpe-Head(-)/Body(+) - -
## 7dpe-Head(+)/Body(+) 0.7858 -
## AA-21dpi-sucrose 0.0511 0.0893
## AA-21dpi-blood 0.3882 0.7766
## 21dpe-Head(-)/Body(-) 0.0511 0.2036
## 21dpe-Head(-)/Body(+) 1.0000 0.9000
## 21dpe-Head(+)/Body(+) 0.0199 0.0624
## AA-21dpi-sucrose AA-21dpi-blood 21dpe-Head(-)/Body(-)
## AA-7dpi-blood - - -
## 7dpe-Head(-)/Body(-) - - -
## 7dpe-Head(-)/Body(+) - - -
## 7dpe-Head(+)/Body(+) - - -
## AA-21dpi-sucrose - - -
## AA-21dpi-blood 0.7967 - -
## 21dpe-Head(-)/Body(-) 0.5669 0.7766 -
## 21dpe-Head(-)/Body(+) 0.5669 0.7895 0.5336
## 21dpe-Head(+)/Body(+) 0.5436 0.7895 0.9160
## 21dpe-Head(-)/Body(+)
## AA-7dpi-blood -
## 7dpe-Head(-)/Body(-) -
## 7dpe-Head(-)/Body(+) -
## 7dpe-Head(+)/Body(+) -
## AA-21dpi-sucrose -
## AA-21dpi-blood -
## 21dpe-Head(-)/Body(-) -
## 21dpe-Head(-)/Body(+) -
## 21dpe-Head(+)/Body(+) 0.2493
##
## P value adjustment method: BH1.1.5 PCoA on species level
#### 7 dpe #####
GP.M1_species_select_0_7dpi = subset_samples(GP.M1_species_select_0, dpi == "7dpi")
GP.M1_species_select_0_rmZIKV <- subset_taxa(GP.M1_species_select_0_7dpi, !taxa_names(GP.M1_species_select_0_7dpi) == "GIM628_NODE_1_length_10744_cov_806.867254")
GP.M1_species_select_0_rmZIKV## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 39 taxa and 42 samples ]
## sample_data() Sample Data: [ 42 samples by 19 sample variables ]
## tax_table() Taxonomy Table: [ 39 taxa by 7 taxonomic ranks ]sample_variables(GP.M1_species_select_0_rmZIKV)## [1] "mosq" "body"
## [3] "mosq_or" "head"
## [5] "Index" "mosq_species"
## [7] "Treatment" "dpi"
## [9] "group.plaque.qPCR.1000." "group.plaque.qPCR.0."
## [11] "group.plaque.qPCR.1000.body.1000." "group.plaque.qPCR.0.body.0."
## [13] "Info1" "plaque_assay"
## [15] "qPCR_head._ZIKV_WNV" "qPCR_body._ZIKV_WNV"
## [17] "head_ZIKV.1000" "body_ZIKV.1000"
## [19] "hd"GP.ord.GP.M1_species_select_0_rmZIKV <- ordinate(GP.M1_species_select_0_rmZIKV, "PCoA", "bray")
plot_ordination(GP.M1_species_select_0_rmZIKV, GP.ord.GP.M1_species_select_0_rmZIKV, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.1000.body.1000.")+
theme_bw()+
geom_point(size = 4) +
scale_shape_manual(values=rev(shap[2:6])) +
scale_color_manual(values = c(col[2],col[5],col[9]))+
theme(legend.title=element_blank(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of eukaryotic virus in Ae. aegypti at 7 dpe excluding ZIKV") +
stat_ellipse(type = "norm", linetype = 2,level=0.6,aes(group = Treatment))
set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_species_select_0_rmZIKV, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_species_select_0_rmZIKV))
# Adonis test
adonis(GP.M1_species_select_bray ~ group.plaque.qPCR.1000.body.1000., data = GP.M1_species_select_sampledf)##
## Call:
## adonis(formula = GP.M1_species_select_bray ~ group.plaque.qPCR.1000.body.1000., data = GP.M1_species_select_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## group.plaque.qPCR.1000.body.1000. 4 3.2397 0.80992 2.5814 0.21818 0.001
## Residuals 37 11.6087 0.31375 0.78182
## Total 41 14.8484 1.00000
##
## group.plaque.qPCR.1000.body.1000. ***
## Residuals
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### 21 dpe #####
GP.M1_species_select_0_21dpi = subset_samples(GP.M1_species_select_0, dpi == "21dpi")
GP.M1_species_select_0_rmZIKV <- subset_taxa(GP.M1_species_select_0_21dpi, !taxa_names(GP.M1_species_select_0_21dpi) == "GIM628_NODE_1_length_10744_cov_806.867254")
GP.M1_species_select_0_rmZIKV## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 39 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 19 sample variables ]
## tax_table() Taxonomy Table: [ 39 taxa by 7 taxonomic ranks ]GP.ord.GP.M1_species_select_0_rmZIKV <- ordinate(GP.M1_species_select_0_rmZIKV, "PCoA", "bray")
plot_ordination(GP.M1_species_select_0_rmZIKV, GP.ord.GP.M1_species_select_0_rmZIKV, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.1000.body.1000.")+
theme_bw()+
geom_point(size = 4) +
scale_shape_manual(values=c(shap[1:4],shap[6],shap[5])) +
scale_color_manual(values = c(col[2],col[5],col[9]))+
theme(legend.title=element_blank(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of eukaryotic virus in Ae. aegypti at 21 dpe excluding ZIKV")
set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_species_select_0_rmZIKV, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_species_select_0_rmZIKV))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)##
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 2 0.7426 0.37129 0.96082 0.03928 0.473
## Residuals 47 18.1622 0.38643 0.96072
## Total 49 18.9048 1.000001.2 Eukaryotic virome of Culex quinquefasciatus
1.2.1 Subset eukaryotic virome of Culex quinquefasciatus
GP.M1_CQ <- subset_samples(GP.M1, mosq_species == "Cx_quinquefasciatus")
GP.M1_CQ <- subset_taxa(GP.M1_CQ, species != "Kaiowa virus")
GP.M1_CQ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 143 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 143 taxa by 7 taxonomic ranks ]1.2.2 merge on species level
GP.M1_species <- GP.M1_CQ %>%
tax_glom(taxrank = "species")
GP.M1_species## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 60 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 60 taxa by 7 taxonomic ranks ]GP.M1_samselect = prune_samples(sample_sums(GP.M1_species) > 0, GP.M1_species)
GP.M1_samselect ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 60 taxa and 55 samples ]
## sample_data() Sample Data: [ 55 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 60 taxa by 7 taxonomic ranks ]GP.M1_species_select_0 = prune_taxa(taxa_sums(GP.M1_samselect) > 0, GP.M1_samselect)
GP.M1_species_select_0## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 27 taxa and 55 samples ]
## sample_data() Sample Data: [ 55 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 27 taxa by 7 taxonomic ranks ]GP.M1_species_select_0 <- subset_taxa(GP.M1_species_select_0, species != "Aedes aegypti totivirus")
GP.M1_species_select_0## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 26 taxa and 55 samples ]
## sample_data() Sample Data: [ 55 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 26 taxa by 7 taxonomic ranks ]GP.M1_species_select = prune_taxa(taxa_sums(GP.M1_samselect) > 500, GP.M1_samselect)
GP.M1_species_select <- subset_taxa(GP.M1_species_select, species != "Aedes aegypti totivirus")
GP.M1_species_select## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 9 taxa and 55 samples ]
## sample_data() Sample Data: [ 55 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 9 taxa by 7 taxonomic ranks ]1.2.3 Heatmap on viral species level
#### get table of reads for each viral species ####
OTU_species <- merge(otu_table(GP.M1_species_select), tax_table(GP.M1_species_select), by=0)
rownames(OTU_species) <- OTU_species$species
rownames(OTU_species)## [1] "unclassified Chloriridovirus" "West Nile virus"
## [3] "Lysoka partiti-like virus" "unclassified Circoviridae"
## [5] "Guadeloupe Culex tymo-like virus" "Marine protobacilladnavirus 1"
## [7] "Lampyris noctiluca errantivirus 1" "Wenzhou sobemo-like virus 3"
## [9] "Phasi Charoen-like phasivirus"#### order viral species by familiy ####
OTU_species$family <- as.character(OTU_species$family)
OTU_species$family## [1] "Iridoviridae" "Flaviviridae" "Picobirnaviridae"
## [4] "Circoviridae" NA "Bacilladnaviridae"
## [7] "Metaviridae" NA "Phenuiviridae"OTU_species$family[is.na(OTU_species$family)] <- "unclassified"
OTU_species_m<- melt(OTU_species)## Using Row.names, superkingdom, phylum, class, order, family, genus, species as id variablesOTU_species_m_agg <- aggregate(value~family+species, FUN=sum, data=OTU_species_m)
OTU_species_m_agg1 <- aggregate(value~family, FUN=sum, data=OTU_species_m)
family_ord <- OTU_species_m_agg1[order(-OTU_species_m_agg1$value),]$family
family_ord <- c(family_ord[-1],"unclassified")
family_ord## [1] "Flaviviridae" "Phenuiviridae" "Metaviridae"
## [4] "Circoviridae" "Picobirnaviridae" "Iridoviridae"
## [7] "Bacilladnaviridae" "unclassified"OTU_species_m_agg$index <- 1
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Phenuiviridae")] <- 2
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Metaviridae")] <- 3
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Circoviridae")] <- 4
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Picobirnaviridae")] <- 5
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Iridoviridae")] <- 6
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "Bacilladnaviridae")] <- 7
OTU_species_m_agg$index[which(OTU_species_m_agg$family == "unclassified")] <- 8
OTU_species_m_agg <- OTU_species_m_agg[order(OTU_species_m_agg$index, -OTU_species_m_agg$value),]
OTU_species_m_agg## family species value index
## 9 Flaviviridae West Nile virus 93203 1
## 5 Phenuiviridae Phasi Charoen-like phasivirus 8886 2
## 2 Metaviridae Lampyris noctiluca errantivirus 1 1874 3
## 7 Circoviridae unclassified Circoviridae 1597 4
## 3 Picobirnaviridae Lysoka partiti-like virus 699 5
## 6 Iridoviridae unclassified Chloriridovirus 575 6
## 4 Bacilladnaviridae Marine protobacilladnavirus 1 503 7
## 1 unclassified Guadeloupe Culex tymo-like virus 4168168 8
## 8 unclassified Wenzhou sobemo-like virus 3 1066803 8species_ord <- as.character(OTU_species_m_agg$species)
species_ord## [1] "West Nile virus" "Phasi Charoen-like phasivirus"
## [3] "Lampyris noctiluca errantivirus 1" "unclassified Circoviridae"
## [5] "Lysoka partiti-like virus" "unclassified Chloriridovirus"
## [7] "Marine protobacilladnavirus 1" "Guadeloupe Culex tymo-like virus"
## [9] "Wenzhou sobemo-like virus 3"which(species_ord == "West Nile virus")## [1] 1species_ord_1 <- c("West Nile virus", species_ord[-1])
species_ord_1## [1] "West Nile virus" "Phasi Charoen-like phasivirus"
## [3] "Lampyris noctiluca errantivirus 1" "unclassified Circoviridae"
## [5] "Lysoka partiti-like virus" "unclassified Chloriridovirus"
## [7] "Marine protobacilladnavirus 1" "Guadeloupe Culex tymo-like virus"
## [9] "Wenzhou sobemo-like virus 3"#### transform the table - reads for each viral species #####
rm <- c("Row.names", "superkingdom", "phylum", "class", "order", "family", "genus", "species")
OTU_species <- OTU_species[,!colnames(OTU_species) %in% rm]
#log
OTU_species_trans <- log10(OTU_species)
is.na(OTU_species_trans) <- sapply(OTU_species_trans, is.infinite)
OTU_species_trans[is.na(OTU_species_trans)]<-0
df.original <- OTU_species
df <- OTU_species_trans
#### use meta data for top annotation #####
meta = read.csv("~/Desktop/mosq_Guad_infection/GIM_metadata_new.csv", header=TRUE, row.names=1, sep=",")
meta_s = meta[!meta$mosq_species == "",]
meta_CQ <- meta_s[colnames(df),]
#### order sample by WNV reads within each group & column annotation ####
meta_CQz <- merge(meta_CQ, as.data.frame(t(df)),by=0)
meta_CQz$qPCR_body._ZIKV_WNV <- log10(meta_CQz$qPCR_body._ZIKV_WNV)
which(colnames(meta_CQz)=="qPCR_body._ZIKV_WNV")## [1] 17is.na(meta_CQz[,17])<-sapply(meta_CQz[,17], is.infinite)
meta_CQz[,17][is.na(meta_CQz[,17])]<-0
meta_CQz[,17]## [1] 4.209086 5.472123 3.452553 3.541704 4.729788 3.716671 4.211761 0.000000
## [9] 2.972203 4.242566 3.451479 0.000000 3.004321 0.000000 2.136721 4.681802
## [17] 0.000000 3.691258 3.618048 4.242566 4.335879 0.000000 0.000000 4.521177
## [25] 4.574945 0.000000 4.350326 5.005627 0.000000 4.113609 3.787744 0.000000
## [33] 0.000000 4.256092 2.775246 4.221336 0.000000 0.000000 0.000000 0.000000
## [41] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## [49] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000meta_CQz$plaque_assay## [1] "7dpi-WNV(-)" "7dpi-WNV(-)" "7dpi-WNV(-)" "7dpi-WNV(+)"
## [5] "7dpi-WNV(+)" "7dpi-WNV(-)" "7dpi-WNV(-)" "7dpi-WNV(+)"
## [9] "7dpi-WNV(-)" "7dpi-WNV(-)" "7dpi-WNV(-)" "7dpi-WNV(-)"
## [13] "7dpi-WNV(-)" "7dpi-WNV(-)" "7dpi-WNV(-)" "7dpi-WNV(-)"
## [17] "14dpi-WNV(+)" "14dpi-WNV(-)" "14dpi-WNV(-)" "14dpi-WNV(-)"
## [21] "14dpi-WNV(-)" "14dpi-WNV(-)" "14dpi-WNV(-)" "14dpi-WNV(+)"
## [25] "14dpi-WNV(+)" "14dpi-WNV(-)" "14dpi-WNV(-)" "14dpi-WNV(+)"
## [29] "14dpi-WNV(-)" "14dpi-WNV(-)" "14dpi-WNV(+)" "14dpi-WNV(-)"
## [33] "14dpi-WNV(-)" "14dpi-WNV(-)" "14dpi-WNV(-)" "14dpi-WNV(-)"
## [37] "CQ-7dpi-blood" "CQ-7dpi-blood" "CQ-7dpi-blood" "CQ-7dpi-blood"
## [41] "CQ-7dpi-blood" "CQ-14dpi-blood" "CQ-14dpi-blood" "CQ-14dpi-blood"
## [45] "CQ-14dpi-blood" "CQ-14dpi-blood" "CQ-7dpi-sucrose" "CQ-7dpi-sucrose"
## [49] "CQ-7dpi-sucrose" "CQ-7dpi-sucrose" "CQ-14dpi-sucrose" "CQ-14dpi-sucrose"
## [53] "CQ-14dpi-sucrose" "CQ-14dpi-sucrose" "CQ-14dpi-sucrose"meta_CQz$plaque_assay <- factor(meta_CQz$plaque_assay,
levels = c("CQ-7dpi-sucrose", "CQ-7dpi-blood", "7dpi-WNV(-)", "7dpi-WNV(+)", "CQ-14dpi-sucrose", "CQ-14dpi-blood", "14dpi-WNV(-)", "14dpi-WNV(+)"))
meta_CQz_order1 <- meta_CQz[meta_CQz$dpi== "7dpi",][with(meta_CQz[meta_CQz$dpi== "7dpi",],
order(plaque_assay, -qPCR_body._ZIKV_WNV, -`West Nile virus`, -`Guadeloupe Culex tymo-like virus`,
-`Wenzhou sobemo-like virus 3`, -`Phasi Charoen-like phasivirus`)), ]
meta_CQz_order2 <- meta_CQz[meta_CQz$dpi== "14dpi",][with(meta_CQz[meta_CQz$dpi== "14dpi",],
order(plaque_assay, -qPCR_body._ZIKV_WNV, -`West Nile virus`, -`Guadeloupe Culex tymo-like virus`,
-`Wenzhou sobemo-like virus 3`, -`Phasi Charoen-like phasivirus`)), ]
meta_CQz_order <- rbind(meta_CQz_order1,meta_CQz_order2)
meta_CQz$Row.names## [1] "GIM382" "GIM383" "GIM384" "GIM385" "GIM387" "GIM388" "GIM390" "GIM392"
## [9] "GIM393" "GIM394" "GIM395" "GIM396" "GIM397" "GIM398" "GIM399" "GIM400"
## [17] "GIM421" "GIM422" "GIM423" "GIM424" "GIM425" "GIM426" "GIM427" "GIM428"
## [25] "GIM429" "GIM430" "GIM431" "GIM432" "GIM433" "GIM434" "GIM435" "GIM436"
## [33] "GIM437" "GIM438" "GIM439" "GIM440" "GIM451" "GIM452" "GIM453" "GIM454"
## [41] "GIM455" "GIM471" "GIM472" "GIM473" "GIM475" "GIM476" "GIM492" "GIM493"
## [49] "GIM494" "GIM495" "GIM525" "GIM527" "GIM528" "GIM529" "GIM530"meta_CQz_order$Row.names## [1] "GIM495" "GIM493" "GIM494" "GIM492" "GIM452" "GIM451" "GIM453" "GIM454"
## [9] "GIM455" "GIM383" "GIM400" "GIM394" "GIM390" "GIM382" "GIM388" "GIM384"
## [17] "GIM395" "GIM397" "GIM393" "GIM399" "GIM398" "GIM396" "GIM387" "GIM385"
## [25] "GIM392" "GIM525" "GIM527" "GIM529" "GIM530" "GIM528" "GIM471" "GIM475"
## [33] "GIM473" "GIM472" "GIM476" "GIM431" "GIM425" "GIM438" "GIM424" "GIM440"
## [41] "GIM434" "GIM422" "GIM423" "GIM439" "GIM426" "GIM436" "GIM433" "GIM430"
## [49] "GIM427" "GIM437" "GIM432" "GIM429" "GIM428" "GIM435" "GIM421"df <- df[,meta_CQz_order$Row.names]
col_ha = HeatmapAnnotation(qPCR = anno_lines(meta_CQz_order[,17], gp = gpar(col = c("darkgrey")),
height = unit(1.5, "cm"),ylim = c(-0.5,6), axis_param = list(at = c(0,2,4,6)),
add_points = TRUE, pt_gp = gpar(col = c("darkgrey")), pch = c(17)))
meta_CQz_order$group_n1 ## NULLmeta_CQz_order$plaque_assay## [1] CQ-7dpi-sucrose CQ-7dpi-sucrose CQ-7dpi-sucrose CQ-7dpi-sucrose
## [5] CQ-7dpi-blood CQ-7dpi-blood CQ-7dpi-blood CQ-7dpi-blood
## [9] CQ-7dpi-blood 7dpi-WNV(-) 7dpi-WNV(-) 7dpi-WNV(-)
## [13] 7dpi-WNV(-) 7dpi-WNV(-) 7dpi-WNV(-) 7dpi-WNV(-)
## [17] 7dpi-WNV(-) 7dpi-WNV(-) 7dpi-WNV(-) 7dpi-WNV(-)
## [21] 7dpi-WNV(-) 7dpi-WNV(-) 7dpi-WNV(+) 7dpi-WNV(+)
## [25] 7dpi-WNV(+) CQ-14dpi-sucrose CQ-14dpi-sucrose CQ-14dpi-sucrose
## [29] CQ-14dpi-sucrose CQ-14dpi-sucrose CQ-14dpi-blood CQ-14dpi-blood
## [33] CQ-14dpi-blood CQ-14dpi-blood CQ-14dpi-blood 14dpi-WNV(-)
## [37] 14dpi-WNV(-) 14dpi-WNV(-) 14dpi-WNV(-) 14dpi-WNV(-)
## [41] 14dpi-WNV(-) 14dpi-WNV(-) 14dpi-WNV(-) 14dpi-WNV(-)
## [45] 14dpi-WNV(-) 14dpi-WNV(-) 14dpi-WNV(-) 14dpi-WNV(-)
## [49] 14dpi-WNV(-) 14dpi-WNV(-) 14dpi-WNV(+) 14dpi-WNV(+)
## [53] 14dpi-WNV(+) 14dpi-WNV(+) 14dpi-WNV(+)
## 8 Levels: CQ-7dpi-sucrose CQ-7dpi-blood 7dpi-WNV(-) ... 14dpi-WNV(+)#### column split ####
df_1 <- df
df_1$WNV <- "others"
df_1$WNV[which(rownames(df_1)=="West Nile virus")] <- "aWNV"#### draw figure #####
heatmap_col_tp2 <- colorRampPalette(brewer.pal(9, "YlGnBu"), space = "rgb")(100)
Heatmap(as.matrix(df),
name = "log10 reads number", #title of legend
column_title = "Samples", row_title = "viral species",
row_names_gp = gpar(fontsize = 10),
column_names_gp = gpar(fontsize = 6),
column_title_gp = gpar(fontsize = 12),
cluster_rows = F,
cluster_columns= F,
row_order = species_ord_1,
column_labels = meta_CQz_order$mosq,
top_annotation = col_ha,
col = heatmap_col_tp2,
column_split = meta_CQz_order$plaque_assay,
row_split = df_1$WNV,
border = TRUE,
rect_gp = gpar(col = "white", lwd = 0.5)) 
2 qRT-PCR analysis
copies = read.table("~/Desktop/mosq_Guad_infection/qRT-PCR/qRTPCR_copies.csv", header=TRUE, dec=".", sep=",")
colnames(copies)## [1] "mosq_No" "mosq_species"
## [3] "Treatment" "dpi"
## [5] "part" "Body"
## [7] "Head" "PCLPV_copies"
## [9] "GMV_copies" "ATV_copies"
## [11] "AANV_copies" "ZIKV_copies"
## [13] "GCTLV_copies" "WSLV3_copies"
## [15] "WNV_copies" "group.plaque.qPCR.1000."
## [17] "group.plaque.qPCR.0." "group.plaque.qPCR.1000.body.1000."
## [19] "group.plaque.qPCR.0.body.0." "Info1"
## [21] "plaque_assay" "head_ZIKV.1000"
## [23] "body_ZIKV.1000"#### log transform ####
copies_log <- copies
copies_log <- cbind(copies[,1:7], log10(copies[,8:15]),copies[,16:23])
is.na(copies_log) <- sapply(copies_log, is.infinite)
copies_log[is.na(copies_log)] <- 0
colnames(copies_log)## [1] "mosq_No" "mosq_species"
## [3] "Treatment" "dpi"
## [5] "part" "Body"
## [7] "Head" "PCLPV_copies"
## [9] "GMV_copies" "ATV_copies"
## [11] "AANV_copies" "ZIKV_copies"
## [13] "GCTLV_copies" "WSLV3_copies"
## [15] "WNV_copies" "group.plaque.qPCR.1000."
## [17] "group.plaque.qPCR.0." "group.plaque.qPCR.1000.body.1000."
## [19] "group.plaque.qPCR.0.body.0." "Info1"
## [21] "plaque_assay" "head_ZIKV.1000"
## [23] "body_ZIKV.1000"########### all core viruses #########
copies_log$plaque_assay <- gsub("CQ-14dpi-","CQ-",copies_log$plaque_assay,fixed = T)
copies_log$plaque_assay <- gsub("CQ-7dpi-","CQ-",copies_log$plaque_assay,fixed = T)
copies_log$plaque_assay <- gsub("14dpi-","CQ-",copies_log$plaque_assay,fixed = T)
copies_log$plaque_assay <- gsub("7dpi-","CQ-",copies_log$plaque_assay,fixed = T)
copies_log$plaque_assay <- factor(copies_log$plaque_assay, levels = c("AA-sucrose", "AA-blood", "AA-ZIKV(-)", "AA-ZIKV(+)",
"CQ-sucrose", "CQ-blood", "CQ-WNV(-)", "CQ-WNV(+)"))
levels(copies_log$plaque_assay)## [1] "AA-sucrose" "AA-blood" "AA-ZIKV(-)" "AA-ZIKV(+)" "CQ-sucrose"
## [6] "CQ-blood" "CQ-WNV(-)" "CQ-WNV(+)"copies_log$PCLPV_copies[copies_log$mosq_species == "Cx_quinquefasciatus"] <- NA
copies_log$GMV_copies[copies_log$mosq_species == "Cx_quinquefasciatus"] <- NA
copies_log$ATV_copies[copies_log$mosq_species == "Cx_quinquefasciatus"] <- NA
copies_log$AANV_copies[copies_log$mosq_species == "Cx_quinquefasciatus"] <- NA
copies_log$GCTLV_copies[copies_log$mosq_species == "Ae_aegypti"] <- NA
copies_log$WSLV3_copies[copies_log$mosq_species == "Ae_aegypti"] <- NA
copies_log_m <- melt(copies_log)## Using mosq_No, mosq_species, Treatment, dpi, part, Body, Head, group.plaque.qPCR.1000., group.plaque.qPCR.0., group.plaque.qPCR.1000.body.1000., group.plaque.qPCR.0.body.0., Info1, plaque_assay, head_ZIKV.1000, body_ZIKV.1000 as id variablescopies_log_m_sub <- copies_log_m[copies_log_m$variable %in% c("PCLPV_copies", "GMV_copies", "ATV_copies","AANV_copies","GCTLV_copies", "WSLV3_copies"),]
copies_log_m_sub <- na.omit(copies_log_m_sub)
copies_log_m_sub$variable <- factor(copies_log_m_sub$variable, levels=c("GMV_copies", "PCLPV_copies","ATV_copies","AANV_copies","GCTLV_copies", "WSLV3_copies"))
levels(copies_log_m_sub$variable)## [1] "GMV_copies" "PCLPV_copies" "ATV_copies" "AANV_copies" "GCTLV_copies"
## [6] "WSLV3_copies"copies_log_m_sub$virus[copies_log_m_sub$variable == "GMV_copies"] <- "Guadeloupe mosquito virus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "PCLPV_copies"] <- "Phasi Charoen−like phasivirus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "ATV_copies"] <- "Aedes aegypti totivirus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "AANV_copies"] <- "Aedes anphevirus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "GCTLV_copies"] <- "Guadeloupe Culex tymo−like virus"
copies_log_m_sub$virus[copies_log_m_sub$variable == "WSLV3_copies"] <- "Wenzhou sobemo−like virus 3"
copies_log_m_sub$virus <- factor(copies_log_m_sub$virus, levels=c("Guadeloupe mosquito virus", "Phasi Charoen−like phasivirus","Aedes aegypti totivirus",
"Aedes anphevirus","Guadeloupe Culex tymo−like virus", "Wenzhou sobemo−like virus 3"))
copies_log_m_sub$dpi <- gsub("dpi", "dpe", copies_log_m_sub$dpi)
copies_log_m_sub$dpi<- factor(copies_log_m_sub$dpi, levels = c("7dpe", "21dpe", "14dpe"))
copies_log_m_sub_rmATV21 <- copies_log_m_sub[!(copies_log_m_sub$dpi=="21dpe"&copies_log_m_sub$variable=="ATV_copies"), ]
copies_log_m_sub_rmGCTLV7 <- copies_log_m_sub_rmATV21[!(copies_log_m_sub_rmATV21$dpi=="7dpe"&copies_log_m_sub_rmATV21$variable=="GCTLV_copies"), ]
copies_log_m_sub_rmWSLV37 <- copies_log_m_sub_rmGCTLV7[!(copies_log_m_sub_rmGCTLV7$dpi=="7dpe"&copies_log_m_sub_rmGCTLV7$variable=="WSLV3_copies"), ]
col = plasma(8,alpha = 1, begin=0.9, end = 0, direction = 1)ggplot(copies_log_m_sub_rmWSLV37, aes(x=plaque_assay, y=value, col=plaque_assay)) +
geom_boxplot(na.rm = TRUE,outlier.shape = NA) + geom_jitter(shape=19, position=position_jitter(0.25), size = 2.2, alpha=0.7)+
theme_bw()+
facet_wrap(~ virus+part, scales="free", ncol = 6)+
scale_color_manual(values = col) +
theme(legend.position = "right",
axis.ticks.x = element_blank(),
axis.title.y = element_text(size=10), axis.title.x = element_blank(),
axis.text.x = element_blank(),
strip.text.x = element_text(size = 12, colour = "black"),
plot.title = element_text(hjust = 0.5, size = 12))+
ggtitle("Viral genome copies determined by qRT-PCR")+
ylab("log10 copies")+
scale_y_continuous(limits = c(-0.01,9), breaks=seq(0, 9, by = 1))
copies_log_m_sub_rmATV7 <- copies_log_m_sub[!(copies_log_m_sub$dpi=="7dpe"&copies_log_m_sub$variable=="ATV_copies"), ]
copies_log_m_sub_rmGCTLV21 <- copies_log_m_sub_rmATV7[!(copies_log_m_sub_rmATV7$dpi=="14dpe"&copies_log_m_sub_rmATV7$variable=="GCTLV_copies"), ]
copies_log_m_sub_rmWSLV321 <- copies_log_m_sub_rmGCTLV21[!(copies_log_m_sub_rmGCTLV21$dpi=="14dpe"&copies_log_m_sub_rmGCTLV21$variable=="WSLV3_copies"), ]
copies_log_m_sub_ATV <- copies_log_m_sub[copies_log_m_sub$variable=="ATV_copies", ]
copies_log_m_sub_ATV$dpi <- "7dpe+21dpe"
copies_log_m_sub_GCTLV <- copies_log_m_sub[copies_log_m_sub$variable=="GCTLV_copies", ]
copies_log_m_sub_GCTLV$dpi <- "7dpe+14dpe"
copies_log_m_sub_WSLV3 <- copies_log_m_sub[copies_log_m_sub$variable=="WSLV3_copies", ]
copies_log_m_sub_WSLV3$dpi <- "7dpe+14dpe"
copies_log_m_sub_new <- rbind(copies_log_m_sub_rmWSLV321,copies_log_m_sub_ATV)
copies_log_m_sub_new <- rbind(copies_log_m_sub_new,copies_log_m_sub_GCTLV)
copies_log_m_sub_new <- rbind(copies_log_m_sub_new,copies_log_m_sub_WSLV3)
ggplot(copies_log_m_sub_new, aes(x=plaque_assay, y=value, col=plaque_assay)) +
geom_boxplot(na.rm = TRUE,outlier.shape = NA) + geom_jitter(shape=19, position=position_jitter(0.25), size = 2.2, alpha=0.7)+
theme_bw()+
facet_wrap(~ virus+dpi+part, scales="free", ncol = 4)+
scale_color_manual(values = col) +
theme(legend.position = "right",
axis.ticks.x = element_blank(),
axis.title.y = element_text(size=10), axis.title.x = element_blank(),
axis.text.x = element_blank(),
strip.text.x = element_text(size = 12, colour = "black"),
plot.title = element_text(hjust = 0.5, size = 12))+
ggtitle("Viral genome copies determined by qRT-PCR")+
ylab("log10 copies")+
scale_y_continuous(limits = c(-0.01,9), breaks=seq(0, 9, by = 1))
########## pairwise.wilcox.test #######
#### GMV Guadeloupe mosquito virus
pairwise.wilcox.test(copies_log$GMV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$GMV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"]
##
## AA-sucrose AA-blood
## AA-blood 0.04491 -
## AA-ZIKV(-) 1.2e-05 0.00059
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$GMV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$GMV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"]
##
## AA-sucrose AA-blood
## AA-blood 0.0720 -
## AA-ZIKV(-) 0.0011 0.0011
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$GMV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$GMV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"]
##
## AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood 0.0507 - -
## AA-ZIKV(-) 0.6847 0.0022 -
## AA-ZIKV(+) 0.2667 0.0225 0.2667
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$GMV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$GMV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"]
##
## AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood 0.0507 - -
## AA-ZIKV(-) 0.0764 0.0031 -
## AA-ZIKV(+) 0.3714 0.0225 0.9538
##
## P value adjustment method: BH#### PCLPV Phasi Charoen−like phasivirus
pairwise.wilcox.test(copies_log$PCLPV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: copies_log$PCLPV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"]
##
## AA-sucrose AA-blood
## AA-blood 1.00 -
## AA-ZIKV(-) 0.44 0.44
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$PCLPV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$PCLPV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"]
##
## AA-sucrose AA-blood
## AA-blood 0.98 -
## AA-ZIKV(-) 0.98 0.98
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$PCLPV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"],
p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: copies_log$PCLPV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"]
##
## AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood 0.84 - -
## AA-ZIKV(-) 0.35 0.32 -
## AA-ZIKV(+) 0.35 0.32 0.46
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$PCLPV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$PCLPV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"]
##
## AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood 0.94 - -
## AA-ZIKV(-) 0.94 0.94 -
## AA-ZIKV(+) 0.94 0.94 0.94
##
## P value adjustment method: BH#### ATV Aedes aegypti totivirus
pairwise.wilcox.test(copies_log$ATV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$ATV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"]
##
## AA-sucrose AA-blood
## AA-blood 0.11 -
## AA-ZIKV(-) 0.03 0.52
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$ATV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$ATV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"]
##
## AA-sucrose AA-blood
## AA-blood 0.100 -
## AA-ZIKV(-) 0.061 1.000
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$ATV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$ATV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"]
##
## AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood 0.93 - -
## AA-ZIKV(-) 0.93 0.93 -
## AA-ZIKV(+) 0.93 0.93 0.93
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$ATV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$ATV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"]
##
## AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood 0.359 - -
## AA-ZIKV(-) 0.366 0.039 -
## AA-ZIKV(+) 0.824 0.216 0.366
##
## P value adjustment method: BH#### AANV Aedes anphevirus
pairwise.wilcox.test(copies_log$AANV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$AANV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "body"]
##
## AA-sucrose AA-blood
## AA-blood 0.84 -
## AA-ZIKV(-) 0.84 0.84
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$AANV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$AANV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "7dpi" & copies_log$part == "head"]
##
## AA-sucrose AA-blood
## AA-blood 1 -
## AA-ZIKV(-) 1 1
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$AANV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$AANV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "body"]
##
## AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood 1.00 - -
## AA-ZIKV(-) 0.58 0.68 -
## AA-ZIKV(+) 0.58 0.58 0.58
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$AANV_copies[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Ae_aegypti" & copies_log$dpi=="21dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$AANV_copies[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Ae_aegypti" & copies_log$dpi == "21dpi" & copies_log$part == "head"]
##
## AA-sucrose AA-blood AA-ZIKV(-)
## AA-blood 1.0 - -
## AA-ZIKV(-) 0.9 0.9 -
## AA-ZIKV(+) 0.9 0.9 0.9
##
## P value adjustment method: BH#### GCTLV Guadeloupe Culex tymo−like virus
pairwise.wilcox.test(copies_log$GCTLV_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: copies_log$GCTLV_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "body"]
##
## CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood 0.62 - -
## CQ-WNV(-) 0.62 1.00 -
## CQ-WNV(+) 1.00 0.82 0.62
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$GCTLV_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$GCTLV_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "head"]
##
## CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood 0.553 - -
## CQ-WNV(-) 0.264 0.097 -
## CQ-WNV(+) 0.264 0.048 0.553
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$GCTLV_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: copies_log$GCTLV_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "body"]
##
## CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood 1.00 - -
## CQ-WNV(-) 0.59 0.59 -
## CQ-WNV(+) 1.00 1.00 0.59
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$GCTLV_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="head"],
p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: copies_log$GCTLV_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "head"]
##
## CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood 0.421 - -
## CQ-WNV(-) 0.011 0.011 -
## CQ-WNV(+) 0.016 0.048 0.202
##
## P value adjustment method: BH#### WSLV3 Wenzhou sobemo−like virus 3
pairwise.wilcox.test(copies_log$WSLV3_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$WSLV3_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "body"]
##
## CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood 0.18 - -
## CQ-WNV(-) 0.68 0.18 -
## CQ-WNV(+) 0.51 0.15 0.18
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$WSLV3_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="7dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$WSLV3_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "7dpi" & copies_log$part == "head"]
##
## CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood 0.90 - -
## CQ-WNV(-) 0.90 0.90 -
## CQ-WNV(+) 0.85 0.85 0.85
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$WSLV3_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="body"],
copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="body"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$WSLV3_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "body"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "body"]
##
## CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood 0.242 - -
## CQ-WNV(-) 0.032 0.794 -
## CQ-WNV(+) 0.797 0.149 0.032
##
## P value adjustment method: BHpairwise.wilcox.test(copies_log$WSLV3_copies[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="head"],
copies_log$plaque_assay[copies_log$mosq_species=="Cx_quinquefasciatus" & copies_log$dpi=="14dpi" & copies_log$part=="head"],
p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: copies_log$WSLV3_copies[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "head"] and copies_log$plaque_assay[copies_log$mosq_species == "Cx_quinquefasciatus" & copies_log$dpi == "14dpi" & copies_log$part == "head"]
##
## CQ-sucrose CQ-blood CQ-WNV(-)
## CQ-blood 0.447 - -
## CQ-WNV(-) 0.064 0.230 -
## CQ-WNV(+) 0.797 0.447 0.064
##
## P value adjustment method: BH3 16S rRNA analysis
#### load phyloseq project ####
GIMBac_decont <- readRDS("~/Desktop/mosq_Guad_infection/16s_rRNA/GIM_Bacteriome_decont_SLV132.rds")
GIMBac_decont## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 503 taxa and 158 samples ]
## sample_data() Sample Data: [ 158 samples by 20 sample variables ]
## tax_table() Taxonomy Table: [ 503 taxa by 6 taxonomic ranks ]col_6_Ae <-c("#8D634A", "#CCC591", "#84A0A5", "#972D15", "#C27D38", "#1D769C")
col_6_Cq <-c("#7A7358", "#85A789", "#C6CDF7","#29211F", "#376138","#7294D4")3.1 Bacteriome of Aedes aegypti
GIMBac.Aa <- subset_samples(GIMBac_decont, mosq_species == "Ae_aegypti")
GIMBac.Aa <- subset_taxa(GIMBac.Aa, taxa_sums(GIMBac.Aa) > 0 )
GIMBac.Aa## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 381 taxa and 94 samples ]
## sample_data() Sample Data: [ 94 samples by 20 sample variables ]
## tax_table() Taxonomy Table: [ 381 taxa by 6 taxonomic ranks ]3.1.1 Genera proportion per sample
#### collape on Genus ####
GIMBac.Aa.Genus <- GIMBac.Aa %>% #raw counts
tax_glom(taxrank = "Genus")
GIMBac.Aa.Genus.per <- GIMBac.Aa.Genus %>%
transform_sample_counts(function(x) {100*(x/sum(x))} ) %>%
psmelt()
GIMBac.Aa.Genus.per$group.plaque.qPCR.1000.body.1000.[which(is.na(GIMBac.Aa.Genus.per$group.plaque.qPCR.1000.body.1000.))] <- "NC"
levels(GIMBac.Aa.Genus.per$plaque_assay)## NULL#### 把一些低丰度的anno改成"others" ####
AaBac_agg <- aggregate(Abundance ~ mosq + plaque_assay + Genus, FUN=sum, data=GIMBac.Aa.Genus.per)
AaBac_wide <- reshape2::dcast(AaBac_agg, mosq + plaque_assay ~ Genus, value.var = 'Abundance')
AaBac_wide <- AaBac_wide[order(-AaBac_wide$g_Asaia),]
sample_order <- AaBac_wide$mosq
AaBac_name_10 <- colnames(AaBac_wide[,3:ncol(AaBac_wide)][,colSums(AaBac_wide[,3:ncol(AaBac_wide)])<10])
nc <- which(colnames(GIMBac.Aa.Genus.per) == "Genus")
nc## [1] 29GIMBac.Aa.Genus.per[GIMBac.Aa.Genus.per$Genus %in% AaBac_name_10, nc] <- "Others"
GIMBac.Aa.Genus.per$Genus <- as.factor(GIMBac.Aa.Genus.per$Genus)
levels(GIMBac.Aa.Genus.per$Genus)## [1] "g_Acidovorax" "g_Acinetobacter"
## [3] "g_Asaia" "g_Bacillus"
## [5] "g_Brevundimonas" "g_Chryseobacterium"
## [7] "g_Cupriavidus" "g_Klebsiella"
## [9] "g_Leucobacter" "g_Mycobacterium"
## [11] "g_Novosphingobium" "g_Pseudomonas"
## [13] "g_Salmonella" "g_Sphingobacterium"
## [15] "g_Sphingomonas" "Others"
## [17] "uc_f_Enterobacteriaceae"###--- repeat ---###
AaBac_agg <- aggregate(Abundance ~ mosq + plaque_assay + Genus, FUN=sum, data=GIMBac.Aa.Genus.per)
AaBac_wide <- reshape2::dcast(AaBac_agg, mosq + plaque_assay ~ Genus, value.var = 'Abundance')
as.data.frame(sort(colSums(AaBac_wide[,3:ncol(AaBac_wide)]) , decreasing=T))## sort(colSums(AaBac_wide[, 3:ncol(AaBac_wide)]), decreasing = T)
## g_Asaia 6866.75066
## g_Chryseobacterium 784.51212
## g_Pseudomonas 509.97242
## g_Novosphingobium 371.14382
## g_Acidovorax 250.81813
## g_Klebsiella 206.27917
## g_Acinetobacter 102.45281
## Others 79.11557
## g_Salmonella 45.36808
## uc_f_Enterobacteriaceae 44.29459
## g_Cupriavidus 30.49819
## g_Mycobacterium 29.13343
## g_Sphingomonas 19.51334
## g_Bacillus 17.45504
## g_Leucobacter 16.68823
## g_Sphingobacterium 14.56847
## g_Brevundimonas 11.43594order <- rev(rownames(as.data.frame(sort(colSums(AaBac_wide[,3:ncol(AaBac_wide)]) , decreasing=T))))
which(order == "Others" | order == "uc_f_Enterobacteriaceae")## [1] 8 10order.1 <- c("Others", "uc_f_Enterobacteriaceae", order[-which(order == "Others" | order == "uc_f_Enterobacteriaceae")])
order.1## [1] "Others" "uc_f_Enterobacteriaceae"
## [3] "g_Brevundimonas" "g_Sphingobacterium"
## [5] "g_Leucobacter" "g_Bacillus"
## [7] "g_Sphingomonas" "g_Mycobacterium"
## [9] "g_Cupriavidus" "g_Salmonella"
## [11] "g_Acinetobacter" "g_Klebsiella"
## [13] "g_Acidovorax" "g_Novosphingobium"
## [15] "g_Pseudomonas" "g_Chryseobacterium"
## [17] "g_Asaia"head(AaBac_agg)## mosq plaque_assay Genus Abundance
## 1 AaB-21dpi-1 AA-21dpi-blood g_Acidovorax 0
## 2 AaB-21dpi-2 AA-21dpi-blood g_Acidovorax 0
## 3 AaB-21dpi-3 AA-21dpi-blood g_Acidovorax 0
## 4 AaB-21dpi-4 AA-21dpi-blood g_Acidovorax 0
## 5 AaB-21dpi-5 AA-21dpi-blood g_Acidovorax 0
## 6 AaS-21dpi-1 AA-21dpi-sucrose g_Acidovorax 0AaBac_agg$Genus <- factor(AaBac_agg$Genus, levels = order.1)
AaBac_agg$plaque_assay <- factor(AaBac_agg$plaque_assay,
levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV(-)","AA-21dpi-sucrose", "AA-21dpi-blood", "AA-21dpi-ZIKV(-)", "AA-21dpi-ZIKV(+)"))
AaBac_agg$mosq <- factor(AaBac_agg$mosq, levels = sample_order)
col <- c("#CBD588", "#5F7FC7", "orange","#DA5724", "#508578", "#CD9BCD",
"#AD6F3B", "#673770","#D14285", "#652926", "#C84248",
"#8569D5", "#5E738F","#D1A33D", "#8A7C64", "#599861","grey")
ggplot(AaBac_agg, aes(x = mosq, y = Abundance, fill = Genus)) +
facet_grid(~ plaque_assay,scales="free",space = "free") + #### Lay out panels in a grid
geom_bar(stat = "identity")+
scale_fill_manual(values = rev(col)) +
theme_bw()+
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5,hjust = 1),
axis.title.y = element_text(size = 12,vjust = 0.5),
plot.margin = unit(c(1, 1, 1, 1), "cm"),
strip.text.x = element_text(size =8, colour = "black"),
plot.title = element_text(hjust = 0.5))
3.1.2 Alpha diversity on genus level
#### rarefy counts ####
sort(rowSums(otu_table(GIMBac.Aa)))## GIM-320 GIM-318 GIM-319 GIM-317 GIM-549 GIM-30 GIM-257 GIM-259 GIM-580 GIM-586
## 787 1445 1653 3635 4538 4978 9910 11428 16039 19805
## GIM-34 GIM-543 GIM-37 GIM-553 GIM-31 GIM-628 GIM-614 GIM-609 GIM-350 GIM-99
## 20068 20530 21983 28947 33256 36063 42868 44769 45373 47031
## GIM-97 GIM-572 GIM-98 GIM-92 GIM-575 GIM-616 GIM-86 GIM-21 GIM-94 GIM-100
## 47309 48020 48458 48841 50160 50288 50498 51337 51395 51467
## GIM-627 GIM-618 GIM-617 GIM-36 GIM-260 GIM-611 GIM-26 GIM-95 GIM-24 GIM-32
## 51962 52368 52446 52954 53403 53839 54581 55520 55666 55936
## GIM-83 GIM-82 GIM-87 GIM-574 GIM-290 GIM-84 GIM-29 GIM-582 GIM-585 GIM-81
## 56574 56680 57223 57323 57427 57896 58005 58291 58441 59844
## GIM-38 GIM-347 GIM-546 GIM-39 GIM-621 GIM-288 GIM-35 GIM-610 GIM-89 GIM-625
## 59877 60555 60803 61113 61342 61460 62137 62399 62488 62557
## GIM-287 GIM-23 GIM-578 GIM-615 GIM-33 GIM-612 GIM-622 GIM-623 GIM-573 GIM-620
## 63010 63549 63662 63813 64146 65068 65498 65642 67201 67355
## GIM-284 GIM-583 GIM-576 GIM-577 GIM-258 GIM-22 GIM-608 GIM-581 GIM-96 GIM-28
## 67642 67652 68252 68459 68582 68584 68685 68986 69076 70193
## GIM-349 GIM-88 GIM-626 GIM-25 GIM-554 GIM-619 GIM-93 GIM-27 GIM-85 GIM-91
## 70253 71987 72692 72716 73162 73723 74380 74467 74633 76212
## GIM-613 GIM-90 GIM-40 GIM-624
## 76585 82853 83026 90098set.seed(100)
GIMBac.Aa_rare = rarefy_even_depth(GIMBac.Aa, sample.size = 11428, replace=FALSE, trimOTUs = F)## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`## ...## 7 samples removedbecause they contained fewer reads than `sample.size`.## Up to first five removed samples are:## GIM-318GIM-30GIM-319GIM-320GIM-549## ...GIMBac.Aa_rm = prune_samples(sample_sums(GIMBac.Aa)<11428, GIMBac.Aa)
GIMBac.Aa_rm## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 381 taxa and 7 samples ]
## sample_data() Sample Data: [ 7 samples by 20 sample variables ]
## tax_table() Taxonomy Table: [ 381 taxa by 6 taxonomic ranks ]GIMBac.Aa_rareall <- merge_phyloseq(GIMBac.Aa_rare, GIMBac.Aa_rm)
sort(rowSums(otu_table(GIMBac.Aa_rareall)))## GIM-320 GIM-318 GIM-319 GIM-317 GIM-549 GIM-30 GIM-257 GIM-21 GIM-29 GIM-36
## 787 1445 1653 3635 4538 4978 9910 11428 11428 11428
## GIM-575 GIM-585 GIM-87 GIM-94 GIM-609 GIM-617 GIM-625 GIM-543 GIM-22 GIM-37
## 11428 11428 11428 11428 11428 11428 11428 11428 11428 11428
## GIM-576 GIM-586 GIM-88 GIM-95 GIM-610 GIM-618 GIM-626 GIM-284 GIM-23 GIM-38
## 11428 11428 11428 11428 11428 11428 11428 11428 11428 11428
## GIM-577 GIM-81 GIM-89 GIM-96 GIM-611 GIM-619 GIM-627 GIM-287 GIM-24 GIM-31
## 11428 11428 11428 11428 11428 11428 11428 11428 11428 11428
## GIM-39 GIM-578 GIM-82 GIM-97 GIM-612 GIM-620 GIM-628 GIM-288 GIM-25 GIM-32
## 11428 11428 11428 11428 11428 11428 11428 11428 11428 11428
## GIM-40 GIM-580 GIM-83 GIM-90 GIM-98 GIM-613 GIM-621 GIM-290 GIM-347 GIM-26
## 11428 11428 11428 11428 11428 11428 11428 11428 11428 11428
## GIM-33 GIM-572 GIM-581 GIM-84 GIM-91 GIM-99 GIM-614 GIM-622 GIM-258 GIM-349
## 11428 11428 11428 11428 11428 11428 11428 11428 11428 11428
## GIM-27 GIM-34 GIM-573 GIM-582 GIM-85 GIM-92 GIM-100 GIM-615 GIM-623 GIM-259
## 11428 11428 11428 11428 11428 11428 11428 11428 11428 11428
## GIM-546 GIM-350 GIM-28 GIM-35 GIM-574 GIM-583 GIM-86 GIM-93 GIM-608 GIM-616
## 11428 11428 11428 11428 11428 11428 11428 11428 11428 11428
## GIM-624 GIM-260 GIM-553 GIM-554
## 11428 11428 11428 11428#### collape on Genus ####
GIMBac.Aa_rare.Genus <- GIMBac.Aa_rareall %>% #raw counts
tax_glom(taxrank = "Genus")
#### alpha diversity plot ####
plot_richness <- function (physeq, x = "samples", color = NULL, shape = NULL,
title = NULL, scales = "free_y", nrow = 1, shsi = NULL, measures = NULL,
sortby = NULL)
{
erDF = estimate_richness(physeq, split = TRUE, measures = measures)
measures = colnames(erDF)
ses = colnames(erDF)[grep("^se\\.", colnames(erDF))]
measures = measures[!measures %in% ses]
if (!is.null(sample_data(physeq, errorIfNULL = FALSE))) {
DF <- data.frame(erDF, sample_data(physeq))
}
else {
DF <- data.frame(erDF)
}
if (!"samples" %in% colnames(DF)) {
DF$samples <- sample_names(physeq)
}
if (!is.null(x)) {
if (x %in% c("sample", "samples", "sample_names", "sample.names")) {
x <- "samples"
}
}
else {
x <- "samples"
}
mdf = reshape2::melt(DF, measure.vars = measures)
mdf$se <- NA_integer_
if (length(ses) > 0) {
selabs = ses
names(selabs) <- substr(selabs, 4, 100)
substr(names(selabs), 1, 1) <- toupper(substr(names(selabs),
1, 1))
mdf$wse <- sapply(as.character(mdf$variable), function(i,
selabs) {
selabs[i]
}, selabs)
for (i in 1:nrow(mdf)) {
if (!is.na(mdf[i, "wse"])) {
mdf[i, "se"] <- mdf[i, (mdf[i, "wse"])]
}
}
mdf <- mdf[, -which(colnames(mdf) %in% c(selabs, "wse"))]
}
if (!is.null(measures)) {
if (any(measures %in% as.character(mdf$variable))) {
mdf <- mdf[as.character(mdf$variable) %in% measures,
]
}
else {
warning("Argument to `measures` not supported. All alpha-diversity measures (should be) included in plot.")
}
}
if (!is.null(shsi)) {
warning("shsi no longer supported option in plot_richness. Please use `measures` instead")
}
if (!is.null(sortby)) {
if (!all(sortby %in% levels(mdf$variable))) {
warning("`sortby` argument not among `measures`. Ignored.")
}
if (!is.discrete(mdf[, x])) {
warning("`sortby` argument provided, but `x` not a discrete variable. `sortby` is ignored.")
}
if (all(sortby %in% levels(mdf$variable)) & is.discrete(mdf[,
x])) {
wh.sortby = which(mdf$variable %in% sortby)
mdf[, x] <- factor(mdf[, x], levels = names(sort(tapply(X = mdf[wh.sortby,
"value"], INDEX = mdf[wh.sortby, x], mean, na.rm = TRUE,
simplify = TRUE))))
}
}
richness_map = aes_string(x = x, y = "value", colour = color,
shape = shape)
p = ggplot(mdf, richness_map) + geom_boxplot(na.rm = TRUE,outlier.shape = NA)+geom_jitter(shape=19, position=position_jitter(0.25), size = 2.2, alpha=0.7)
# if (any(!is.na(mdf[, "se"]))) {
# p = p + geom_errorbar(aes(ymax = value + se, ymin = value -
# se), width = 0.1)
# }
p = p + theme(axis.text.x = element_text(angle = -90, vjust = 0.5,
hjust = 0))
p = p + ylab("Alpha Diversity Measure")
p = p + facet_wrap(~variable, nrow = nrow, scales = scales)
if (!is.null(title)) {
p <- p + ggtitle(title)
}
return(p)
}
sample_data(GIMBac.Aa_rare.Genus)$Info1 <- factor(sample_data(GIMBac.Aa_rare.Genus)$Info1, levels = c("AA-7dpi-sucrose", "AA-7dpi-blood", "AA-7dpi-ZIKV",
"AA-21dpi-sucrose","AA-21dpi-blood", "AA-21dpi-ZIKV"))
sample_data(GIMBac.Aa_rare.Genus)$Treatment <- factor(sample_data(GIMBac.Aa_rare.Genus)$Treatment, levels = c("AA-sucrose", "AA-blood", "AA-ZIKV"))
sample_data(GIMBac.Aa_rare.Genus)$dpi <- factor(sample_data(GIMBac.Aa_rare.Genus)$dpi, levels = c("7dpi", "21dpi"))
GIMBac.Aa_rare.Genus.7d = subset_samples(GIMBac.Aa_rare.Genus, dpi == "7dpi")
GIMBac.Aa_rare.Genus.21d = subset_samples(GIMBac.Aa_rare.Genus, dpi == "21dpi")
plot_richness(GIMBac.Aa_rare.Genus, x="Info1",measures=c("Shannon"),col="Treatment", shape = "dpi") +
scale_color_manual(values = rev(col_6_Ae[c(4:6)])) +
theme_bw() +
ylim(-0.001,1.9)+
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 12, colour = "black", angle = 90, hjust = 1, vjust=0.5),
strip.text.x = element_text(size = 12, colour = "black"),
plot.title = element_text(hjust = 0.5, size = 12)) +
ggtitle("Alpha diversity of bacteriome in Ae. aegypti on genus level")+
ylab("Alpha Diversity Value")
erich <- estimate_richness(GIMBac.Aa_rare.Genus, measures = c("Shannon"))
pairwise.wilcox.test(erich$Shannon,sample_data(GIMBac.Aa_rare.Genus)$Info1, p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: erich$Shannon and sample_data(GIMBac.Aa_rare.Genus)$Info1
##
## AA-7dpi-sucrose AA-7dpi-blood AA-7dpi-ZIKV AA-21dpi-sucrose
## AA-7dpi-blood 0.0433 - - -
## AA-7dpi-ZIKV 0.2523 0.0075 - -
## AA-21dpi-sucrose 0.0433 0.5867 0.0433 -
## AA-21dpi-blood 0.0170 0.0298 3.0e-05 0.0170
## AA-21dpi-ZIKV 0.0028 0.9455 2.3e-08 0.2916
## AA-21dpi-blood
## AA-7dpi-blood -
## AA-7dpi-ZIKV -
## AA-21dpi-sucrose -
## AA-21dpi-blood -
## AA-21dpi-ZIKV 0.0028
##
## P value adjustment method: BH3.1.3 Beta diversity on genus level
GIMBac.Aa_rare.GenusS <- subset_samples(GIMBac.Aa_rare.Genus, group.plaque.qPCR.1000.body.1000. != "")
sample_data(GIMBac.Aa_rare.GenusS)$group.plaque.qPCR.1000.body.1000. <- factor(sample_data(GIMBac.Aa_rare.GenusS)$group.plaque.qPCR.1000.body.1000.,
levels =c("AA-7dpi-ZIKV(-/+/+)", "AA-7dpi-ZIKV(-/-/+)", "AA-7dpi-ZIKV(-/-/-)",
"AA-7dpi-blood", "AA-7dpi-sucrose", "AA-21dpi-ZIKV(+/+/+)",
"AA-21dpi-ZIKV(-/+/+)", "AA-21dpi-ZIKV(-/-/+)", "AA-21dpi-ZIKV(-/-/-)",
"AA-21dpi-blood", "AA-21dpi-sucrose"))
GIMBac.Aa_rare.Genus.7dS = subset_samples(GIMBac.Aa_rare.GenusS, dpi == "7dpi")
GIMBac.Aa_rare.Genus.21dS = subset_samples(GIMBac.Aa_rare.GenusS, dpi == "21dpi")#### 7dpi #####
levels(sample_data(GIMBac.Aa_rare.Genus.7dS)$group.plaque.qPCR.1000.body.1000.)## [1] "AA-7dpi-ZIKV(-/+/+)" "AA-7dpi-ZIKV(-/-/+)" "AA-7dpi-ZIKV(-/-/-)"
## [4] "AA-7dpi-blood" "AA-7dpi-sucrose"levels(sample_data(GIMBac.Aa_rare.Genus.7dS)$Treatment)## [1] "AA-sucrose" "AA-blood" "AA-ZIKV"GP.ord.GIMBac.Aa_rare.Genus.7d <- ordinate(GIMBac.Aa_rare.Genus.7dS, "PCoA", "bray")
plot_ordination(GIMBac.Aa_rare.Genus.7dS, GP.ord.GIMBac.Aa_rare.Genus.7d, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.1000.body.1000.")+
theme_bw()+
geom_point(size = 4) +
scale_color_manual(values = col_6_Ae[c(4:6)])+
scale_shape_manual(values=rev(shap[2:6])) +
theme(legend.title=element_blank(),
# legend.text= element_text(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of bacteriome in Ae. aegypti at 7 dpe") +
stat_ellipse(type = "norm", linetype = 2,level=0.6,aes(group = Treatment))
set.seed(1)
# Calculate bray curtis distance matrix
GIMBac.Aa_rare.Genus.7dS_bray <- phyloseq::distance(GIMBac.Aa_rare.Genus.7dS, method = "bray")
# make a data frame from the sample_data
GIMBac.Aa_rare.Genus.7dS_sampledf <- data.frame(sample_data(GIMBac.Aa_rare.Genus.7dS))
# Adonis test
adonis(GIMBac.Aa_rare.Genus.7dS_bray ~ Treatment, data = GIMBac.Aa_rare.Genus.7dS_sampledf)##
## Call:
## adonis(formula = GIMBac.Aa_rare.Genus.7dS_bray ~ Treatment, data = GIMBac.Aa_rare.Genus.7dS_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 2 2.8805 1.44024 10.932 0.35923 0.001 ***
## Residuals 39 5.1381 0.13175 0.64077
## Total 41 8.0186 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### 21dpi #####
levels(sample_data(GIMBac.Aa_rare.Genus.21dS)$group.plaque.qPCR.1000.body.1000.)## [1] "AA-21dpi-ZIKV(+/+/+)" "AA-21dpi-ZIKV(-/+/+)" "AA-21dpi-ZIKV(-/-/+)"
## [4] "AA-21dpi-ZIKV(-/-/-)" "AA-21dpi-blood" "AA-21dpi-sucrose"levels(sample_data(GIMBac.Aa_rare.Genus.21dS)$Treatment)## [1] "AA-sucrose" "AA-blood" "AA-ZIKV"GP.ord.GIMBac.Aa_rare.Genus.21d <- ordinate(GIMBac.Aa_rare.Genus.21dS, "PCoA", "bray")
plot_ordination(GIMBac.Aa_rare.Genus.21dS, GP.ord.GIMBac.Aa_rare.Genus.21d, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.1000.body.1000.")+
theme_bw()+
geom_point(size = 4) +
scale_color_manual(values = col_6_Ae[c(4:6)])+
scale_shape_manual(values= c(shap[1:4],shap[6],shap[5])) +
theme(legend.title=element_blank(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of bacteriome in Ae. aegypti at 21 dpe")
set.seed(1)
# Calculate bray curtis distance matrix
GIMBac.Aa_rare.Genus.21dS_bray <- phyloseq::distance(GIMBac.Aa_rare.Genus.21dS, method = "bray")
# make a data frame from the sample_data
GIMBac.Aa_rare.Genus.21dS_sampledf <- data.frame(sample_data(GIMBac.Aa_rare.Genus.21dS))
# Adonis test
adonis(GIMBac.Aa_rare.Genus.21dS_bray ~ Treatment, data = GIMBac.Aa_rare.Genus.21dS_sampledf)##
## Call:
## adonis(formula = GIMBac.Aa_rare.Genus.21dS_bray ~ Treatment, data = GIMBac.Aa_rare.Genus.21dS_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 2 0.68676 0.34338 8.9531 0.27588 0.002 **
## Residuals 47 1.80259 0.03835 0.72412
## Total 49 2.48935 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.2 Bacteriome of Culex quinquefasciatus
GIMBac.Cq <- subset_samples(GIMBac_decont, mosq_species == "Cx_quinquefasciatus")
GIMBac.Cq <- subset_taxa(GIMBac.Cq, taxa_sums(GIMBac.Cq) > 0 )
GIMBac.Cq## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 154 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 20 sample variables ]
## tax_table() Taxonomy Table: [ 154 taxa by 6 taxonomic ranks ]sample_data(GIMBac.Cq)$group.plaque.qPCR.0.body.0. <- factor(sample_data(GIMBac.Cq)$group.plaque.qPCR.0.body.0.,
levels = c("CQ-7dpi-sucrose","CQ-7dpi-blood","7dpi-WNV(-/-)", "7dpi-WNV(-/+)", "7dpi-WNV(+/)",
"CQ-14dpi-sucrose","CQ-14dpi-blood","14dpi-WNV(-/-)","14dpi-WNV(-/+)","14dpi-WNV(+/)"))
sample_data(GIMBac.Cq)$Info1 <- factor(sample_data(GIMBac.Cq)$Info1,
levels = c("CQ-7dpi-sucrose", "CQ-7dpi-blood", "CQ-7dpi-WNV", "CQ-14dpi-sucrose", "CQ-14dpi-blood", "CQ-14dpi-WNV"))
sample_data(GIMBac.Cq)$Treatment <- factor(sample_data(GIMBac.Cq)$Treatment, levels = c("CQ-sucrose", "CQ-blood", "CQ-WNV"))3.2.1 Genera proportion per sample
#### collape on Genus ####
GIMBac.Cq.Genus <- GIMBac.Cq %>% #raw counts
tax_glom(taxrank = "Genus")
GIMBac.Cq.Genus.per <- GIMBac.Cq.Genus %>%
transform_sample_counts(function(x) {100*(x/sum(x))} ) %>%
psmelt()
colnames(GIMBac.Cq.Genus.per)## [1] "OTU" "Sample"
## [3] "Abundance" "mosq"
## [5] "body" "mosq_or"
## [7] "head" "Index"
## [9] "mosq_species" "Treatment"
## [11] "dpi" "group.plaque.qPCR.1000."
## [13] "group.plaque.qPCR.0." "group.plaque.qPCR.1000.body.1000."
## [15] "group.plaque.qPCR.0.body.0." "Info1"
## [17] "plaque_assay" "qPCR_head._ZIKV_WNV"
## [19] "qPCR_body._ZIKV_WNV" "head_ZIKV.1000"
## [21] "body_ZIKV.1000" "Sample_or_Control"
## [23] "conc" "Kingdom"
## [25] "Phylum" "Class"
## [27] "Order" "Family"
## [29] "Genus"levels(GIMBac.Cq.Genus.per$plaque_assay)## NULLGIMBac.Cq.Genus.per$plaque_assay <- factor(GIMBac.Cq.Genus.per$plaque_assay,
levels = c("CQ-7dpi-sucrose", "CQ-7dpi-blood", "7dpi-WNV(-)", "7dpi-WNV(+)",
"CQ-14dpi-sucrose", "CQ-14dpi-blood", "14dpi-WNV(-)", "14dpi-WNV(+)"))
levels(GIMBac.Cq.Genus.per$plaque_assay)## [1] "CQ-7dpi-sucrose" "CQ-7dpi-blood" "7dpi-WNV(-)" "7dpi-WNV(+)"
## [5] "CQ-14dpi-sucrose" "CQ-14dpi-blood" "14dpi-WNV(-)" "14dpi-WNV(+)"#### 把一些低丰度的anno改成“others” ####
CqBac_agg <- aggregate(Abundance ~ mosq + plaque_assay + Genus, FUN=sum, data=GIMBac.Cq.Genus.per)
CqBac_wide <- dcast(CqBac_agg, mosq + plaque_assay ~ Genus, value.var = 'Abundance')
CqBac_wide_ordered <- CqBac_wide[with(CqBac_wide, order(plaque_assay,-g_Wolbachia,-g_Serratia)), ]
sample_order <- CqBac_wide_ordered$mosq
CqBac_name_10 <- colnames(CqBac_wide[,3:ncol(CqBac_wide)][,colSums(CqBac_wide[,3:ncol(CqBac_wide)])<10])
nc <- which(colnames(GIMBac.Cq.Genus.per) == "Genus")
nc## [1] 29GIMBac.Cq.Genus.per[GIMBac.Cq.Genus.per$Genus %in% CqBac_name_10, nc] <- "Others"
GIMBac.Cq.Genus.per$Genus <- as.factor(GIMBac.Cq.Genus.per$Genus)
levels(GIMBac.Cq.Genus.per$Genus)## [1] "g_Asaia" "g_Delftia" "g_Elizabethkingia"
## [4] "g_Gluconobacter" "g_Pseudomonas" "g_Serratia"
## [7] "g_Stenotrophomonas" "g_Wolbachia" "Others"###--- repeat ---###
CqBac_agg <- aggregate(Abundance ~ mosq + plaque_assay + Genus, FUN=sum, data=GIMBac.Cq.Genus.per)
CqBac_wide <- dcast(CqBac_agg, mosq + plaque_assay ~ Genus, value.var = 'Abundance')
as.data.frame(sort(colSums(CqBac_wide[,3:ncol(CqBac_wide)]) , decreasing=T))## sort(colSums(CqBac_wide[, 3:ncol(CqBac_wide)]), decreasing = T)
## g_Wolbachia 2789.01416
## g_Serratia 1649.54264
## g_Asaia 1244.71360
## g_Elizabethkingia 134.91983
## g_Pseudomonas 71.33818
## g_Gluconobacter 65.51299
## g_Stenotrophomonas 16.96808
## Others 15.88106
## g_Delftia 12.10946order <- rev(rownames(as.data.frame(sort(colSums(CqBac_wide[,3:ncol(CqBac_wide)]) , decreasing=T))))
which(order == "Others")## [1] 2order.1 <- c("Others", order[-which(order == "Others")])
order.1## [1] "Others" "g_Delftia" "g_Stenotrophomonas"
## [4] "g_Gluconobacter" "g_Pseudomonas" "g_Elizabethkingia"
## [7] "g_Asaia" "g_Serratia" "g_Wolbachia"CqBac_agg$mosq <- factor(CqBac_agg$mosq, levels = sample_order)
CqBac_agg$Genus <- factor(CqBac_agg$Genus, levels = order.1)
levels(CqBac_agg$plaque_assay)## [1] "CQ-7dpi-sucrose" "CQ-7dpi-blood" "7dpi-WNV(-)" "7dpi-WNV(+)"
## [5] "CQ-14dpi-sucrose" "CQ-14dpi-blood" "14dpi-WNV(-)" "14dpi-WNV(+)"#### plot Genera proportion bar plot ####
col <- c("#5F7FC7", "#AD6F3B", "#CBD588", "#673770", "orange", "#CD9BCD", "#652926", "#C84248", "grey")
ggplot(CqBac_agg, aes(x = mosq, y = Abundance, fill = Genus)) +
facet_grid(~ plaque_assay,scales="free",space = "free") + #### Lay out panels in a grid
geom_bar(stat = "identity") + #### stat = "identity" : the height of the bar will represent the value in a column of the data frame
scale_fill_manual(values = rev(col)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 8, angle = 90, vjust = 0.5,hjust = 1),
axis.title.y = element_text(size = 12,vjust = 0.5),
plot.margin = unit(c(1, 1, 1, 1), "cm"),
strip.text.x = element_text(size =8, colour = "black"),
plot.title = element_text(hjust = 0.5)) +
guides(fill = guide_legend(reverse = TRUE, keywidth = 1, keyheight = 1))
3.2.2 Proportion of Asaia and Wolbachia
CqBac_agg$dpi <- CqBac_agg$plaque_assay
CqBac_agg$dpi <- gsub("CQ-", "", CqBac_agg$dpi)
CqBac_agg_v <- CqBac_agg %>% separate(dpi, c("dpi", "food"), sep ="-"); head(CqBac_agg_v)## Warning: Expected 2 pieces. Additional pieces discarded in 270 rows [11, 12, 13,
## 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 41, 42, 43, 44, 45, ...].## mosq plaque_assay Genus Abundance dpi food
## 1 CqS-7dpi-1 CQ-7dpi-sucrose g_Asaia 15.4421010 7dpi sucrose
## 2 CqS-7dpi-2 CQ-7dpi-sucrose g_Asaia 6.1206321 7dpi sucrose
## 3 CqS-7dpi-3 CQ-7dpi-sucrose g_Asaia 1.2262223 7dpi sucrose
## 4 CqS-7dpi-4 CQ-7dpi-sucrose g_Asaia 3.4048961 7dpi sucrose
## 5 CqS-7dpi-5 CQ-7dpi-sucrose g_Asaia 1.9477277 7dpi sucrose
## 6 CqB-7dpi-1 CQ-7dpi-blood g_Asaia 0.2025679 7dpi bloodCqBac_agg_v$Genus <- gsub("g_", "", CqBac_agg_v$Genus); head(CqBac_agg_v)## mosq plaque_assay Genus Abundance dpi food
## 1 CqS-7dpi-1 CQ-7dpi-sucrose Asaia 15.4421010 7dpi sucrose
## 2 CqS-7dpi-2 CQ-7dpi-sucrose Asaia 6.1206321 7dpi sucrose
## 3 CqS-7dpi-3 CQ-7dpi-sucrose Asaia 1.2262223 7dpi sucrose
## 4 CqS-7dpi-4 CQ-7dpi-sucrose Asaia 3.4048961 7dpi sucrose
## 5 CqS-7dpi-5 CQ-7dpi-sucrose Asaia 1.9477277 7dpi sucrose
## 6 CqB-7dpi-1 CQ-7dpi-blood Asaia 0.2025679 7dpi bloodCqBac_agg_v <- CqBac_agg_v[CqBac_agg_v$Genus %in% c("Asaia", "Wolbachia"), ]
CqBac_agg_v$dpi <- factor(CqBac_agg_v$dpi, levels= c("7dpi", "14dpi"))
wilcox.test(CqBac_agg_v[CqBac_agg_v$Genus=="Asaia"&CqBac_agg_v$dpi=="7dpi",]$Abundance,
CqBac_agg_v[CqBac_agg_v$Genus=="Asaia"&CqBac_agg_v$dpi=="14dpi",]$Abundance,
p.adjust.method = "BH")##
## Wilcoxon rank sum exact test
##
## data: CqBac_agg_v[CqBac_agg_v$Genus == "Asaia" & CqBac_agg_v$dpi == "7dpi", ]$Abundance and CqBac_agg_v[CqBac_agg_v$Genus == "Asaia" & CqBac_agg_v$dpi == "14dpi", ]$Abundance
## W = 182, p-value = 4.063e-05
## alternative hypothesis: true location shift is not equal to 0wilcox.test(CqBac_agg_v[CqBac_agg_v$Genus=="Wolbachia"&CqBac_agg_v$dpi=="7dpi",]$Abundance,
CqBac_agg_v[CqBac_agg_v$Genus=="Wolbachia"&CqBac_agg_v$dpi=="14dpi",]$Abundance,
p.adjust.method = "BH")##
## Wilcoxon rank sum exact test
##
## data: CqBac_agg_v[CqBac_agg_v$Genus == "Wolbachia" & CqBac_agg_v$dpi == "7dpi", ]$Abundance and CqBac_agg_v[CqBac_agg_v$Genus == "Wolbachia" & CqBac_agg_v$dpi == "14dpi", ]$Abundance
## W = 614, p-value = 0.01487
## alternative hypothesis: true location shift is not equal to 0CqBac_agg_v$Food_sources[CqBac_agg_v$food == "WNV(+)"] <- "WNV(plaque+)"
CqBac_agg_v$Food_sources[CqBac_agg_v$food == "WNV("] <- "WNV(plaque-)"
CqBac_agg_v$Food_sources[CqBac_agg_v$food == "sucrose"] <- "sucrose"
CqBac_agg_v$Food_sources[CqBac_agg_v$food == "blood"] <- "blood"
CqBac_agg_v$Food_sources <- factor(CqBac_agg_v$Food_sources, levels=c("sucrose", "blood", "WNV(plaque-)", "WNV(plaque+)"))
ggplot(CqBac_agg_v, aes(x=Genus, y=Abundance,fill=dpi)) +
geom_boxplot(na.rm = TRUE,outlier.shape = NA)+
facet_wrap(~Food_sources, nrow = 1)+
theme_bw()+
scale_fill_manual(values = col_6_Ae[c(3,1)])+
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 12, hjust = 0.5),
axis.text.y = element_text(size = 12),
axis.title.y = element_text(size = 12,vjust = 0.5),
# plot.margin = unit(c(1, 1, 1, 1), "cm"),
strip.text.x = element_text(size =12, colour = "black"),
plot.title = element_text(hjust = 0.5))+
ylab("Proportion")+
stat_compare_means()
3.2.2 Alpha diversity on genus level
#### rarefy ####
sort(rowSums(otu_table(GIMBac.Cq)))## GIM-434 GIM-495 GIM-530 GIM-527 GIM-440 GIM-400 GIM-452 GIM-382 GIM-396 GIM-392
## 2367 6007 13219 24168 24676 29224 34447 38493 38678 39933
## GIM-491 GIM-493 GIM-435 GIM-472 GIM-529 GIM-436 GIM-388 GIM-525 GIM-394 GIM-387
## 42151 44527 46282 47873 49431 51469 51694 52552 53008 53342
## GIM-473 GIM-528 GIM-391 GIM-433 GIM-475 GIM-454 GIM-381 GIM-430 GIM-437 GIM-492
## 53406 53990 54060 56789 57476 57756 58002 58040 58072 58409
## GIM-397 GIM-425 GIM-432 GIM-453 GIM-399 GIM-384 GIM-393 GIM-383 GIM-389 GIM-427
## 58899 59979 60071 60274 60339 60551 60802 61024 61027 61520
## GIM-424 GIM-471 GIM-428 GIM-431 GIM-421 GIM-451 GIM-385 GIM-390 GIM-422 GIM-438
## 62791 62813 62965 63245 63659 64176 64967 65183 65903 67316
## GIM-476 GIM-426 GIM-429 GIM-395 GIM-423 GIM-398 GIM-455 GIM-494 GIM-386 GIM-439
## 68186 68512 68932 70712 71919 73396 73752 75039 76302 77901set.seed(100)
GIMBac.Cq_rare = rarefy_even_depth(GIMBac.Cq, sample.size = 24168, replace=FALSE, trimOTUs = F)## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`## ...## 3 samples removedbecause they contained fewer reads than `sample.size`.## Up to first five removed samples are:## GIM-495GIM-434GIM-530## ...GIMBac.Cq_rm = prune_samples(sample_sums(GIMBac.Cq)<24168, GIMBac.Cq)
GIMBac.Cq_rm## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 154 taxa and 3 samples ]
## sample_data() Sample Data: [ 3 samples by 20 sample variables ]
## tax_table() Taxonomy Table: [ 154 taxa by 6 taxonomic ranks ]GIMBac.Cq_rareall <- merge_phyloseq(GIMBac.Cq_rare, GIMBac.Cq_rm)
#### collape on Genus ####
GIMBac.Cq_rare.Genus <- GIMBac.Cq_rareall %>% #raw counts
tax_glom(taxrank = "Genus")
#### alpha diversity plot ####
GIMBac.Cq_rare.Genus## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 99 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 20 sample variables ]
## tax_table() Taxonomy Table: [ 99 taxa by 6 taxonomic ranks ]GIMBac.Cq_rare.Genus.7d = subset_samples(GIMBac.Cq_rare.Genus, dpi == "7dpi")
GIMBac.Cq_rare.Genus.14d = subset_samples(GIMBac.Cq_rare.Genus, dpi == "14dpi")
sample_data(GIMBac.Cq_rare.Genus)$dpi <- factor(sample_data(GIMBac.Cq_rare.Genus)$dpi, levels = c("7dpi", "14dpi"))
plot_richness(GIMBac.Cq_rare.Genus, x="Info1",measures=c("Shannon"),col="Treatment", shape = "dpi") +
scale_color_manual(values = rev(col_6_Ae[c(4:6)])) +
theme_bw() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 12, colour = "black", angle = 90, hjust = 1, vjust=0.5),
strip.text.x = element_text(size = 12, colour = "black"),
plot.title = element_text(hjust = 0.5, size = 12)) +
ggtitle("Alpha diversity of bacteriome in Cx. quinquefasciatus on genus level")+
ylab("Alpha Diversity Value")
erich <- estimate_richness(GIMBac.Cq_rare.Genus, measures = c("Shannon"))
pairwise.wilcox.test(erich$Shannon,sample_data(GIMBac.Cq_rare.Genus)$Info1, p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: erich$Shannon and sample_data(GIMBac.Cq_rare.Genus)$Info1
##
## CQ-7dpi-sucrose CQ-7dpi-blood CQ-7dpi-WNV CQ-14dpi-sucrose
## CQ-7dpi-blood 0.32 - - -
## CQ-7dpi-WNV 0.32 0.89 - -
## CQ-14dpi-sucrose 0.86 0.70 0.32 -
## CQ-14dpi-blood 0.58 0.90 0.78 0.32
## CQ-14dpi-WNV 0.97 0.32 0.12 0.78
## CQ-14dpi-blood
## CQ-7dpi-blood -
## CQ-7dpi-WNV -
## CQ-14dpi-sucrose -
## CQ-14dpi-blood -
## CQ-14dpi-WNV 0.14
##
## P value adjustment method: BH3.2.3 Beta diversity on genus level
GIMBac.Cq_rare.Genus## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 99 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 20 sample variables ]
## tax_table() Taxonomy Table: [ 99 taxa by 6 taxonomic ranks ]GIMBac.Cq_rare.Genus.7d## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 99 taxa and 30 samples ]
## sample_data() Sample Data: [ 30 samples by 20 sample variables ]
## tax_table() Taxonomy Table: [ 99 taxa by 6 taxonomic ranks ]GIMBac.Cq_rare.Genus.14d## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 99 taxa and 30 samples ]
## sample_data() Sample Data: [ 30 samples by 20 sample variables ]
## tax_table() Taxonomy Table: [ 99 taxa by 6 taxonomic ranks ]#### 7dpi #####
GP.ord.GIMBac.Cq_rare.Genus.7d <- ordinate(GIMBac.Cq_rare.Genus.7d, "PCoA", "bray")
plot_ordination(GIMBac.Cq_rare.Genus.7d, GP.ord.GIMBac.Cq_rare.Genus.7d, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.0.body.0.")+
theme_bw()+
geom_point(size = 4) +
scale_color_manual(values = col_6_Ae[c(4:6)])+
scale_shape_manual(values=c(shap[1:4],shap[6])) +
theme(legend.title=element_blank(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of bacteriome in Cx. quinquefasciatus at 7dpi")
dev.off()## null device
## 1set.seed(1)
# Calculate bray curtis distance matrix
GIMBac.Cq_rare.Genus.7d_bray <- phyloseq::distance(GIMBac.Cq_rare.Genus.7d, method = "bray")
# make a data frame from the sample_data
GIMBac.Cq_rare.Genus.7d_sampledf <- data.frame(sample_data(GIMBac.Cq_rare.Genus.7d))
# Adonis test
adonis(GIMBac.Cq_rare.Genus.7d_bray ~ Treatment, data = GIMBac.Cq_rare.Genus.7d_sampledf)##
## Call:
## adonis(formula = GIMBac.Cq_rare.Genus.7d_bray ~ Treatment, data = GIMBac.Cq_rare.Genus.7d_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 2 0.5371 0.26853 1.4751 0.0985 0.184
## Residuals 27 4.9153 0.18205 0.9015
## Total 29 5.4524 1.0000#### 14dpi #####
GP.ord.GIMBac.Cq_rare.Genus.14d <- ordinate(GIMBac.Cq_rare.Genus.14d, "PCoA", "bray")
plot_ordination(GIMBac.Cq_rare.Genus.14d, GP.ord.GIMBac.Cq_rare.Genus.14d, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.0.body.0.")+
theme_bw()+
geom_point(size = 4) +
scale_color_manual(values = col_6_Ae[c(4:6)])+
scale_shape_manual(values=c(shap[1:4],shap[6])) +
theme(legend.title=element_blank(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of bacteriome in Cx. quinquefasciatus at 14dpi")
set.seed(1)
# Calculate bray curtis distance matrix
GIMBac.Cq_rare.Genus.14d_bray <- phyloseq::distance(GIMBac.Cq_rare.Genus.14d, method = "bray")
# make a data frame from the sample_data
GIMBac.Cq_rare.Genus.14d_sampledf <- data.frame(sample_data(GIMBac.Cq_rare.Genus.14d))
# Adonis test
adonis(GIMBac.Cq_rare.Genus.14d_bray ~ Treatment, data = GIMBac.Cq_rare.Genus.14d_sampledf)##
## Call:
## adonis(formula = GIMBac.Cq_rare.Genus.14d_bray ~ Treatment, data = GIMBac.Cq_rare.Genus.14d_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 2 0.4523 0.22613 1.1465 0.07828 0.353
## Residuals 27 5.3254 0.19724 0.92172
## Total 29 5.7777 1.000004 DESeq2 determined differential bacteria
#### Aedes ####
df_otu <- as.data.frame(t(GIMBac.Aa_rare.Genus@otu_table))
df_tax <- as.data.frame(GIMBac.Aa_rare.Genus@tax_table)
df_meta <- data.frame(sample_data(GIMBac.Aa_rare.Genus))
df <- merge(df_otu, df_tax, by = "row.names", all=T)
#### Aedes-Diet & Aedes-Infectious ####
AAdiet <- readRDS("~/Desktop/mosq_Guad_infection/RDS/bacteria_DESeq2_AA_treatment.RDS")
AAdiet## [[1]]
## baseMean log2FoldChange lfcSE stat pvalue padj
## SVs5 591.326427 -28.084810 1.525607 -18.408936 1.113902e-75 6.460630e-74
## SVs6 298.034542 -27.219452 1.647940 -16.517262 2.756314e-61 7.993311e-60
## SVs4 1492.902368 -13.216161 1.265541 -10.443093 1.575906e-25 3.046752e-24
## SVs15 35.077201 -8.350071 1.310547 -6.371438 1.872640e-10 2.715327e-09
## SVs12 36.391970 -4.582519 1.002896 -4.569287 4.893853e-06 5.676869e-05
## SVs45 3.197046 -4.904141 1.229969 -3.987206 6.685584e-05 6.462731e-04
## Genus cat
## SVs5 g_Novosphingobium AAbloodvsZIKV
## SVs6 g_Acidovorax AAbloodvsZIKV
## SVs4 g_Chryseobacterium AAbloodvsZIKV
## SVs15 g_Cupriavidus AAbloodvsZIKV
## SVs12 g_Acinetobacter AAbloodvsZIKV
## SVs45 g_Roseococcus AAbloodvsZIKV
##
## [[2]]
## baseMean log2FoldChange lfcSE stat pvalue
## SVs5 495.298867 -23.741984 1.5667411 -15.153738 7.157719e-52
## SVs6 249.635674 -25.546161 1.6928788 -15.090366 1.874008e-51
## SVs9 2277.909223 17.218830 2.1226043 8.112124 4.974231e-16
## SVs4 1254.221186 -5.649115 1.1792208 -4.790549 1.663254e-06
## SVs15 29.380892 -6.098126 1.3126416 -4.645690 3.389420e-06
## SVs78 4.132256 6.416490 1.5232939 4.212247 2.528427e-05
## SVs8 344.495062 2.854758 0.7157120 3.988697 6.643732e-05
## SVs35 9.598872 2.891754 0.8544597 3.384307 7.135821e-04
## SVs102 6.156045 11.388979 3.7627032 3.026808 2.471511e-03
## SVs71 6.554007 5.420686 1.9483073 2.782254 5.398278e-03
## padj Genus cat
## SVs5 2.863088e-50 g_Novosphingobium AAsucrosevsZIKV
## SVs6 3.748015e-50 g_Acidovorax AAsucrosevsZIKV
## SVs9 6.632307e-15 g_Klebsiella AAsucrosevsZIKV
## SVs4 1.663254e-05 g_Chryseobacterium AAsucrosevsZIKV
## SVs15 2.711536e-05 g_Cupriavidus AAsucrosevsZIKV
## SVs78 1.685618e-04 g_Sphingomonas AAsucrosevsZIKV
## SVs8 3.796418e-04 g_Pseudomonas AAsucrosevsZIKV
## SVs35 3.567911e-03 g_Bacillus AAsucrosevsZIKV
## SVs102 1.098449e-02 g_Sphingobacterium AAsucrosevsZIKV
## SVs71 2.159311e-02 g_Brevundimonas AAsucrosevsZIKV
##
## [[3]]
## baseMean log2FoldChange lfcSE stat pvalue padj
## SVs102 22.03441 24.131719 2.8839593 8.367566 5.882061e-17 4.176263e-15
## SVs9 228.68067 8.504409 1.5953683 5.330687 9.784177e-08 3.473383e-06
## SVs8 647.40055 4.005112 0.9860531 4.061761 4.870402e-05 1.152662e-03
## SVs71 19.20429 8.782719 2.2848029 3.843972 1.210586e-04 2.148791e-03
## SVs12 14.66698 5.430173 1.5982411 3.397593 6.798141e-04 9.653360e-03
## Genus cat
## SVs102 g_Sphingobacterium AAsucrosevsblood
## SVs9 g_Klebsiella AAsucrosevsblood
## SVs8 g_Pseudomonas AAsucrosevsblood
## SVs71 g_Brevundimonas AAsucrosevsblood
## SVs12 g_Acinetobacter AAsucrosevsbloodAAinfec <- readRDS("~/Desktop/mosq_Guad_infection/RDS/Bacteria_DESeq2_AA_infectious.RDS")
AAinfec## [[1]]
## baseMean log2FoldChange lfcSE stat pvalue padj
## SVs12 122.6242 7.87385 1.298181 6.065296 1.317112e-09 1.040519e-07
## Genus cat
## SVs12 g_Acinetobacter 21plaqueposvsneggp <- as.data.frame(sort(table(c(AAdiet[[1]]$Genus, AAdiet[[2]]$Genus, AAdiet[[3]]$Genus, AAinfec[[1]]$Genus)), decreasing = T))
gp <- gp[-c(8,10,12),]
order <- c("Chryseobacterium", "Novosphingobium", "Acidovorax", "Cupriavidus", "Roseococcus",
"Acinetobacter", "Klebsiella", "Brevundimonas", "Sphingobacterium")
df_sp <- df[df$Genus %in% as.character(gp$Var1),]
df_sp <- df_sp[, colnames(df_sp) %in% c(colnames(df_otu), "Genus")]
df_spm <- melt(df_sp); head(df_spm)## Using Genus as id variables## Genus variable value
## 1 g_Sphingobacterium GIM-21 0
## 2 g_Acinetobacter GIM-21 17
## 3 g_Cupriavidus GIM-21 42
## 4 g_Chryseobacterium GIM-21 2012
## 5 g_Roseococcus GIM-21 0
## 6 g_Novosphingobium GIM-21 79df_meta$variable <- rownames(df_meta)
df_merge <- merge(df_spm, df_meta, by = "variable", all =T); head(df_merge)## variable Genus value mosq body mosq_or head
## 1 GIM-100 g_Sphingobacterium 0 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 2 GIM-100 g_Acinetobacter 2 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 3 GIM-100 g_Cupriavidus 0 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 4 GIM-100 g_Chryseobacterium 5965 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 5 GIM-100 g_Roseococcus 0 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## 6 GIM-100 g_Novosphingobium 0 AaZ-21dpi-20 GIM-100 AMT-21dpi-20 GIM-80
## Index mosq_species Treatment dpi group.plaque.qPCR.1000.
## 1 AMT-21dpi-B20 Ae_aegypti AA-ZIKV 21dpi AA-21dpi-ZIKV(-/-)
## 2 AMT-21dpi-B20 Ae_aegypti AA-ZIKV 21dpi AA-21dpi-ZIKV(-/-)
## 3 AMT-21dpi-B20 Ae_aegypti AA-ZIKV 21dpi AA-21dpi-ZIKV(-/-)
## 4 AMT-21dpi-B20 Ae_aegypti AA-ZIKV 21dpi AA-21dpi-ZIKV(-/-)
## 5 AMT-21dpi-B20 Ae_aegypti AA-ZIKV 21dpi AA-21dpi-ZIKV(-/-)
## 6 AMT-21dpi-B20 Ae_aegypti AA-ZIKV 21dpi AA-21dpi-ZIKV(-/-)
## group.plaque.qPCR.0. group.plaque.qPCR.1000.body.1000.
## 1 AA-21dpi-ZIKV(-/-) AA-21dpi-ZIKV(-/-/-)
## 2 AA-21dpi-ZIKV(-/-) AA-21dpi-ZIKV(-/-/-)
## 3 AA-21dpi-ZIKV(-/-) AA-21dpi-ZIKV(-/-/-)
## 4 AA-21dpi-ZIKV(-/-) AA-21dpi-ZIKV(-/-/-)
## 5 AA-21dpi-ZIKV(-/-) AA-21dpi-ZIKV(-/-/-)
## 6 AA-21dpi-ZIKV(-/-) AA-21dpi-ZIKV(-/-/-)
## group.plaque.qPCR.0.body.0. Info1 plaque_assay
## 1 AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 2 AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 3 AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 4 AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 5 AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## 6 AA-21dpi-ZIKV(-/-/-) AA-21dpi-ZIKV AA-21dpi-ZIKV(-)
## qPCR_head._ZIKV_WNV qPCR_body._ZIKV_WNV head_ZIKV.1000 body_ZIKV.1000
## 1 0 0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 2 0 0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 3 0 0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 4 0 0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 5 0 0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## 6 0 0 AA-21dpi-ZIKV(-) AA-21dpi-ZIKV(-)
## Sample_or_Control conc
## 1 TRUE 11.0088
## 2 TRUE 11.0088
## 3 TRUE 11.0088
## 4 TRUE 11.0088
## 5 TRUE 11.0088
## 6 TRUE 11.0088df_merge$log10reads <- log10(df_merge$value)
is.na(df_merge$log10reads) <- sapply(df_merge$log10reads, is.infinite)
df_merge$log10reads[is.na(df_merge$log10reads)] <- 0
df_merge$Genus <- gsub("g_", "", df_merge$Genus)
df_merge$Genus <- factor(df_merge$Genus, levels = order); levels(df_merge$Genus)## [1] "Chryseobacterium" "Novosphingobium" "Acidovorax" "Cupriavidus"
## [5] "Roseococcus" "Acinetobacter" "Klebsiella" "Brevundimonas"
## [9] "Sphingobacterium"df_merge_aa <- df_merge
df_merge_aa <- df_merge_aa[,colnames(df_merge_aa) %in% c("Treatment", "Genus", "variable", "log10reads")]
#### Culex ####
df_otu <- as.data.frame(t(GIMBac.Cq_rare.Genus@otu_table))
df_tax <- as.data.frame(GIMBac.Cq_rare.Genus@tax_table)
df_meta <- data.frame(sample_data(GIMBac.Cq_rare.Genus))
df <- merge(df_otu, df_tax, by = "row.names", all=T)
#### Culex-Diet ####
CQdiet <- readRDS("~/Desktop/mosq_Guad_infection/RDS/Bacteria_DESeq2_CQ_treatment.RDS")
CQdiet## [[1]]
## baseMean log2FoldChange lfcSE stat pvalue padj
## SVs3 70296.5 4.432407 1.163301 3.810196 0.0001388567 0.008886832
## Genus cat
## SVs3 g_Serratia CQbloodvsWNV
##
## [[2]]
## baseMean log2FoldChange lfcSE stat pvalue padj
## SVs7 325.366291 30.000000 3.6211361 8.284693 1.184144e-16 9.591566e-15
## SVs85 1.855982 23.153094 3.6296383 6.378898 1.783664e-10 7.223838e-09
## SVs8 305.228196 3.813213 0.7360724 5.180487 2.213078e-07 5.975312e-06
## SVs92 9.928171 8.105716 2.5808287 3.140741 1.685207e-03 3.412545e-02
## Genus cat
## SVs7 g_Elizabethkingia CQsucrosevsWNV
## SVs85 g_Legionella CQsucrosevsWNV
## SVs8 g_Pseudomonas CQsucrosevsWNV
## SVs92 g_Pseudonocardia CQsucrosevsWNV
##
## [[3]]
## baseMean log2FoldChange lfcSE stat pvalue padj
## SVs85 9.093585 23.018909 2.960427 7.775538 7.512755e-15 4.582781e-13
## SVs35 36.037711 6.273055 1.320633 4.750038 2.033788e-06 6.203053e-05
## SVs55 15.306501 7.452474 1.913401 3.894885 9.824549e-05 1.997658e-03
## SVs92 14.322496 7.325056 2.475187 2.959395 3.082438e-03 4.700719e-02
## Genus cat
## SVs85 g_Legionella CQsucrosevsblood
## SVs35 g_Bacillus CQsucrosevsblood
## SVs55 g_Pelomonas CQsucrosevsblood
## SVs92 g_Pseudonocardia CQsucrosevsbloodgp <- as.data.frame(sort(table(c(CQdiet[[1]]$Genus, CQdiet[[2]]$Genus, CQdiet[[3]]$Genus)), decreasing = T))
gp## Var1 Freq
## 1 g_Legionella 2
## 2 g_Pseudonocardia 2
## 3 g_Bacillus 1
## 4 g_Elizabethkingia 1
## 5 g_Pelomonas 1
## 6 g_Pseudomonas 1
## 7 g_Serratia 1order <- c("Serratia", "Pseudomonas", "Bacillus")
df_sp <- df[df$Genus %in% c("g_Serratia", "g_Pseudomonas", "g_Bacillus"),]
df_sp <- df_sp[, colnames(df_sp) %in% c(colnames(df_otu), "Genus")]
df_spm <- melt(df_sp); head(df_spm)## Using Genus as id variables## Genus variable value
## 1 g_Serratia GIM-387 0
## 2 g_Bacillus GIM-387 52
## 3 g_Pseudomonas GIM-387 99
## 4 g_Serratia GIM-394 7252
## 5 g_Bacillus GIM-394 6
## 6 g_Pseudomonas GIM-394 100df_meta$variable <- rownames(df_meta)
df_merge <- merge(df_spm, df_meta, by = "variable", all =T); head(df_merge)## variable Genus value mosq body mosq_or head
## 1 GIM-381 g_Serratia 20620 CqW-7dpi-1 GIM-381 CWV-7dpi-1 GIM-361
## 2 GIM-381 g_Bacillus 5 CqW-7dpi-1 GIM-381 CWV-7dpi-1 GIM-361
## 3 GIM-381 g_Pseudomonas 14 CqW-7dpi-1 GIM-381 CWV-7dpi-1 GIM-361
## 4 GIM-382 g_Serratia 0 CqW-7dpi-2 GIM-382 CWV-7dpi-2 GIM-362
## 5 GIM-382 g_Bacillus 3 CqW-7dpi-2 GIM-382 CWV-7dpi-2 GIM-362
## 6 GIM-382 g_Pseudomonas 110 CqW-7dpi-2 GIM-382 CWV-7dpi-2 GIM-362
## Index mosq_species Treatment dpi group.plaque.qPCR.1000.
## 1 CWV-7dpi-B1 Cx_quinquefasciatus CQ-WNV 7dpi WNV(-)
## 2 CWV-7dpi-B1 Cx_quinquefasciatus CQ-WNV 7dpi WNV(-)
## 3 CWV-7dpi-B1 Cx_quinquefasciatus CQ-WNV 7dpi WNV(-)
## 4 CWV-7dpi-B2 Cx_quinquefasciatus CQ-WNV 7dpi WNV(-)
## 5 CWV-7dpi-B2 Cx_quinquefasciatus CQ-WNV 7dpi WNV(-)
## 6 CWV-7dpi-B2 Cx_quinquefasciatus CQ-WNV 7dpi WNV(-)
## group.plaque.qPCR.0. group.plaque.qPCR.1000.body.1000.
## 1
## 2
## 3
## 4
## 5
## 6
## group.plaque.qPCR.0.body.0. Info1 plaque_assay qPCR_head._ZIKV_WNV
## 1 7dpi-WNV(-/+) CQ-7dpi-WNV 7dpi-WNV(-) 0
## 2 7dpi-WNV(-/+) CQ-7dpi-WNV 7dpi-WNV(-) 0
## 3 7dpi-WNV(-/+) CQ-7dpi-WNV 7dpi-WNV(-) 0
## 4 7dpi-WNV(-/+) CQ-7dpi-WNV 7dpi-WNV(-) 0
## 5 7dpi-WNV(-/+) CQ-7dpi-WNV 7dpi-WNV(-) 0
## 6 7dpi-WNV(-/+) CQ-7dpi-WNV 7dpi-WNV(-) 0
## qPCR_body._ZIKV_WNV head_ZIKV.1000 body_ZIKV.1000 Sample_or_Control conc
## 1 2399 TRUE 21.9990
## 2 2399 TRUE 21.9990
## 3 2399 TRUE 21.9990
## 4 16184 TRUE 10.7685
## 5 16184 TRUE 10.7685
## 6 16184 TRUE 10.7685df_merge$log10reads <- log10(df_merge$value)
is.na(df_merge$log10reads) <- sapply(df_merge$log10reads, is.infinite)
df_merge$log10reads[is.na(df_merge$log10reads)] <- 0
df_merge$Genus <- gsub("g_", "", df_merge$Genus)
df_merge$Genus <- factor(df_merge$Genus, levels = order); levels(df_merge$Genus)## [1] "Serratia" "Pseudomonas" "Bacillus"df_merge_cq <- df_merge
df_merge_cq <- df_merge_cq[,colnames(df_merge_cq) %in% c("Treatment", "Genus", "variable", "log10reads")]
#### Aedes & Culex ####
df_merge_all <- rbind(df_merge_aa, df_merge_cq)
df_merge_all$Treatment## [1] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [7] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [13] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [19] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [25] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [31] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [37] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [43] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [49] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [55] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [61] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [67] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [73] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [79] AA-blood AA-blood AA-blood AA-ZIKV AA-ZIKV AA-ZIKV
## [85] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [91] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [97] AA-blood AA-blood AA-blood AA-ZIKV AA-ZIKV AA-ZIKV
## [103] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [109] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [115] AA-ZIKV AA-ZIKV AA-ZIKV AA-blood AA-blood AA-blood
## [121] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [127] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [133] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [139] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [145] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [151] AA-ZIKV AA-ZIKV AA-ZIKV AA-blood AA-blood AA-blood
## [157] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [163] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [169] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [175] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [181] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [187] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [193] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [199] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [205] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV AA-ZIKV AA-ZIKV
## [211] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [217] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [223] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV AA-ZIKV AA-ZIKV
## [229] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [235] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [241] AA-ZIKV AA-ZIKV AA-ZIKV AA-sucrose AA-sucrose AA-sucrose
## [247] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [253] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [259] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV AA-ZIKV AA-ZIKV
## [265] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [271] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [277] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV AA-ZIKV AA-ZIKV
## [283] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [289] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [295] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [301] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [307] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [313] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [319] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [325] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [331] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [337] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [343] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [349] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [355] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [361] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [367] AA-sucrose AA-sucrose AA-sucrose AA-ZIKV AA-ZIKV AA-ZIKV
## [373] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [379] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [385] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [391] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [397] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [403] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [409] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [415] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [421] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [427] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [433] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [439] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [445] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [451] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [457] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [463] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [469] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [475] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [481] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [487] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [493] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [499] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [505] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [511] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [517] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [523] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [529] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [535] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [541] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [547] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [553] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [559] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [565] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [571] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [577] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [583] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [589] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [595] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [601] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [607] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [613] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [619] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [625] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [631] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [637] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [643] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [649] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [655] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [661] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [667] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [673] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [679] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [685] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [691] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [697] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [703] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [709] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [715] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [721] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [727] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [733] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [739] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [745] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [751] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [757] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [763] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [769] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [775] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [781] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [787] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [793] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [799] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [805] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [811] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [817] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [823] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [829] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [835] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [841] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [847] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [853] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [859] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [865] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [871] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [877] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [883] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [889] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [895] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [901] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [907] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [913] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [919] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [925] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [931] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [937] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [943] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [949] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [955] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [961] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [967] CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood
## [973] CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood
## [979] CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood
## [985] CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood
## [991] CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood
## [997] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [1003] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [1009] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [1015] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [1021] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## Levels: AA-sucrose AA-blood AA-ZIKV CQ-sucrose CQ-blood CQ-WNVdf_merge_all$Treatment <- factor(df_merge_all$Treatment , levels = c("AA-sucrose", "AA-blood", "AA-ZIKV", "CQ-sucrose", "CQ-blood", "CQ-WNV"))#### plot ####
ggplot(df_merge_all, aes(x = Treatment, y = log10reads, fill = Treatment)) +
facet_wrap(~Genus,ncol = 3, scales="free_x") +
theme_bw() +
geom_boxplot(na.rm = TRUE,outlier.shape = NA, alpha = 0.7)+
geom_jitter(shape=19, position=position_jitter(0.2), size = 2,alpha = 0.5)+
scale_fill_manual(values = c(rev(col_6_Ae[c(4:6)]), rev(col_6_Ae[c(4:6)])))+
theme(axis.title.y = element_text(size = 14, vjust=0.5),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
legend.position = "right",
strip.text.x =element_text(size = 14, hjust = 0.5)) +
ylab("log10 rarefied counts")
5 Phageome analysis
GuadPhage <- readRDS("~/Desktop/mosq_Guad_infection/RDS/GIM_phagome_raw.RDS")5.1 HEATMAP of phageome
#### load phyloseq project ####
Phage_host = readRDS("~/Desktop/mosq_Guad_infection/RDS/GIM_phage_host.RDS")
df <- as.data.frame(otu_table(Phage_host))
df_tax <- as.data.frame(tax_table(Phage_host))
df_meta <- data.frame(sample_data(Phage_host))
cat <- c("AA-sucrose", "AA-blood", "AA-ZIKV", "CQ-sucrose", "CQ-blood", "CQ-WNV")
#### mat1 log10 reads per treatment #####
Phage_host_treatment <- Phage_host %>%
tax_glom(taxrank = "Treatment")
Phage_host_treatment_merge <- merge(as.data.frame(Phage_host_treatment@otu_table), as.data.frame(Phage_host_treatment@tax_table), by="row.names")
rownames(Phage_host_treatment_merge) <- Phage_host_treatment_merge$Treatment
Phage_host_treatment_merge <- Phage_host_treatment_merge[, !colnames(Phage_host_treatment_merge) %in% c("Row.names", "mosq_species", "Treatment","dpi", "Info1")]
mat1 <- log10(Phage_host_treatment_merge)
is.na(mat1) <- sapply(mat1, is.infinite)
mat1[is.na(mat1)]<-0
mat1[1:5,1:5]## GIM319_NODE_68_length_3999_cov_6.501275
## AA-blood 3.611723
## AA-sucrose 3.623559
## CQ-WNV 0.000000
## CQ-blood 0.000000
## CQ-sucrose 0.000000
## GIM424_NODE_178_length_2803_cov_5.668379
## AA-blood 2.209515
## AA-sucrose 0.000000
## CQ-WNV 3.046105
## CQ-blood 0.000000
## CQ-sucrose 0.000000
## GIM573_NODE_150_length_3346_cov_8.985317
## AA-blood 3.247237
## AA-sucrose 3.092370
## CQ-WNV 0.000000
## CQ-blood 0.000000
## CQ-sucrose 0.000000
## GIM95_NODE_22_length_3970_cov_5.511174
## AA-blood 3.411620
## AA-sucrose 3.232742
## CQ-WNV 0.000000
## CQ-blood 0.000000
## CQ-sucrose 0.000000
## GIM400_NODE_1_length_17961_cov_34.782767
## AA-blood 2.460898
## AA-sucrose 0.000000
## CQ-WNV 6.087870
## CQ-blood 5.796710
## CQ-sucrose 4.686609mat1 <- mat1[cat,]
rownames(mat1)## [1] "AA-sucrose" "AA-blood" "AA-ZIKV" "CQ-sucrose" "CQ-blood"
## [6] "CQ-WNV"#### mat2 presence% per treatment #####
Phage_host_pre <- Phage_host %>%
tax_glom(taxrank = "Treatment")
Phage_host_pre_merge <- merge(as.data.frame(Phage_host_pre@otu_table), as.data.frame(Phage_host_pre@tax_table), by="row.names")
rownames(Phage_host_pre_merge) <- Phage_host_pre_merge$Treatment
Phage_host_pre_merge <- Phage_host_pre_merge[, !colnames(Phage_host_pre_merge) %in% c("Row.names", "mosq_species", "Treatment","dpi", "Info1")]
mat2 <- Phage_host_pre_merge
mat2 <- mat2[cat,]
rownames(mat2)## [1] "AA-sucrose" "AA-blood" "AA-ZIKV" "CQ-sucrose" "CQ-blood"
## [6] "CQ-WNV"table(df_tax$Treatment)##
## AA-blood AA-sucrose AA-ZIKV CQ-blood CQ-sucrose CQ-WNV
## 10 10 74 10 10 40mat2[1,] <- mat2[1,]/10*100
mat2[2,] <- mat2[2,]/10*100
mat2[3,] <- mat2[3,]/74*100
mat2[4,] <- mat2[4,]/10*100
mat2[5,] <- mat2[5,]/10*100
mat2[6,] <- mat2[6,]/40*100
##abundance log transform
dflog <- log10(df)
is.na(dflog) <- sapply(dflog, is.infinite)
dflog[is.na(dflog)]<-0
#### order rows-sample and columns-contigs ####
##column - phage host
host_or <- data.frame(sort(table(df_meta$host_genus), decreasing=T)); host_or## Var1 Freq
## 1 no result 203
## 2 Wolbachia 23
## 3 Enterobacter 21
## 4 Serratia 8
## 5 Phytobacter 6
## 6 Asaia 2df_meta$host_genus <- factor(df_meta$host_genus, levels = host_or$Var1); levels(df_meta$host_genus)## [1] "no result" "Wolbachia" "Enterobacter" "Serratia" "Phytobacter"
## [6] "Asaia"##column - viral annotation
family_or <- data.frame(sort(table(df_meta$family), decreasing=T)); family_or## Var1 Freq
## 1 no annotation 222
## 2 Siphoviridae 17
## 3 Myoviridae 9
## 4 Podoviridae 7
## 5 Microviridae 6
## 6 Herelleviridae 2df_meta$family <- factor(df_meta$family, levels = family_or$Var1); levels(df_meta$family)## [1] "no annotation" "Siphoviridae" "Myoviridae" "Podoviridae"
## [5] "Microviridae" "Herelleviridae"##column - sample order-length
df_meta$log10length <- log10(df_meta$length)
contigs_ordered <- df_meta[with(df_meta, order(host_genus, -log10length)),]
df_meta$length1k <- df_meta$length/1000
######## draw the map ########
heatmap_col_tp1 <- colorRamp2(c(0:5), brewer.pal(6, "YlOrBr"), space = "RGB") #heatmap colour
heatmap_col_tp2 <- colorRamp2(c(0,20,40,60,80,100), brewer.pal(6, "YlGnBu"), space = "RGB")
heatmap_col_tp3 <- colorRamp2(c(0,20,40,60,80,100), brewer.pal(6, "Greens"), space = "RGB")
#### host block ####
colh = c("#818687FF", "#2D5C21FF", "#2CA0A9FF", "#A1E4BDFF", "#323B34FF", "#C8D656FF")
colhost = HeatmapAnnotation(host = anno_block(gp = gpar(fill = colh)),show_annotation_name = F, annotation_label = F)
#### viral annotation ####
colf = c("#818687FF", "#66C2A5", "#FC8D62", "#8DA0CB", "#FFD92F", "#E5C494")
colfamily = HeatmapAnnotation(condition = df_meta$family,
col = list(condition = c("no annotation"=colf[1], "Siphoviridae"=colf[2], "Myoviridae"=colf[3],
"Podoviridae"=colf[4], "Microviridae"=colf[5], "Herelleviridae"=colf[6]),
show_annotation_name = T, annotation_label = F,
border = T))
#### contigs length ####
collen = HeatmapAnnotation(length = anno_barplot(df_meta$length1k, ylim = c(0,50),
height = unit(2, "cm")))
df_meta$aai_completeness[is.na(df_meta$aai_completeness)] <- 0
complete = HeatmapAnnotation(complete = df_meta$aai_completeness, col = list(complete=heatmap_col_tp3))
#### legend ####
lgd_host = Legend(labels = levels(df_meta$host_genus), title = "Predicted_host", legend_gp = gpar(fill = colh))
lgd_virus = Legend(labels = levels(df_meta$family), title = "Phage_annotation", legend_gp = gpar(fill = colf))
# lgd_reads = Legend(title = "Log10 reads number", col = heatmap_col_tp1, at = c(0:6), labels = c("0", "1", "2", "3", "4", "5", "6"))
# lgd_pro = Legend(title = "Relative presence", col = heatmap_col_tp2, at = c(0:5), labels = c("0", "20", "40", "60", "80", "100"))
pd = packLegend(lgd_host, lgd_virus,
# lgd_reads, lgd_pro,
direction = "vertical")#### plot ####
ht1 <- Heatmap(as.matrix(mat1),
name = "log10 reads number", #title of legend
row_title = "Samples",
col = heatmap_col_tp1,
column_split = df_meta$host_genus,
row_names_gp = gpar(fontsize = 12),
column_names_gp = gpar(fontsize = 0),
column_order = contigs_ordered$name,
row_order = cat,
cluster_rows = F,
cluster_columns= F,
rect_gp = gpar(col = "white", lwd = 0.5))
ht2 <- Heatmap(as.matrix(mat2),
name = "Present proportion", #title of legend
row_title = "Samples",
col = heatmap_col_tp2,
column_split = df_meta$host_genus,
row_names_gp = gpar(fontsize = 12),
column_names_gp = gpar(fontsize = 0),
column_order = contigs_ordered$name,
row_order = cat,
cluster_rows = F,
cluster_columns= F,
rect_gp = gpar(col = "white", lwd = 0.5))
ht_list = colhost %v% colfamily %v% complete %v% collen %v% ht1 %v% ht2
# draw(ht_list, annotation_legend_list = pd)
ht_draw = draw(ht_list, annotation_legend_list = pd)
5.2 Alpha diversity of phageome
GuadPhage_samselect = prune_samples(sample_sums(GuadPhage) > 0 , GuadPhage)
GuadPhage_samselect## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 263 taxa and 150 samples ]
## sample_data() Sample Data: [ 150 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 263 taxa by 17 taxonomic ranks ]factor(sample_data(GuadPhage_samselect)$Treatment)## [1] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [7] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [13] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [19] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [25] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [31] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [37] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [43] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [49] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [55] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [61] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV
## [67] AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-ZIKV AA-sucrose
## [73] AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose AA-sucrose
## [79] AA-sucrose AA-sucrose AA-sucrose AA-blood AA-blood AA-blood
## [85] AA-blood AA-blood AA-blood AA-blood AA-blood AA-blood
## [91] AA-blood CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [97] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [103] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [109] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [115] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [121] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV
## [127] CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-WNV CQ-sucrose
## [133] CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose CQ-sucrose
## [139] CQ-sucrose CQ-sucrose CQ-blood CQ-blood CQ-blood CQ-blood
## [145] CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood CQ-blood
## Levels: AA-blood AA-sucrose AA-ZIKV CQ-blood CQ-sucrose CQ-WNVsample_data(GuadPhage_samselect)$Treatment <- factor(sample_data(GuadPhage_samselect)$Treatment,
levels = c("AA-sucrose","AA-blood", "AA-ZIKV", "CQ-sucrose", "CQ-blood", "CQ-WNV"))
phage <- scan(file = "~/OneDrive/1. Project + Guadeloupe mosq/phage_contigs_complete+20.txt", what="", sep="\n")
my_subset <- subset(otu_table(GuadPhage_samselect), rownames(otu_table(GuadPhage_samselect)) %in% phage)
new_GuadPhage_samselect <- merge_phyloseq(my_subset, tax_table(GuadPhage_samselect), sample_data(GuadPhage_samselect))
plot_richness(new_GuadPhage_samselect, x="Treatment",measures=c("Shannon"),col="Treatment", shape = "mosq_species") +
scale_color_manual(values = rev(c(col[9],col[5],col[2],col[9],col[5],col[2]))) +
theme_bw()+
# facet_wrap(~mosq_species, scales = "free")+
scale_shape_manual(values = c(16, 17))+
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
strip.text.x = element_text(size = 12, colour = "black"),
plot.title = element_text(hjust = 0.5, size = 12),
axis.ticks.x = element_blank()) +
ggtitle("Alpha diversity of phages")+
ylab("Alpha Diversity Value")## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.
erich <- estimate_richness(new_GuadPhage_samselect, measures = c("Shannon"))## Warning in estimate_richness(new_GuadPhage_samselect, measures = c("Shannon")): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.pairwise.wilcox.test(erich$Shannon,sample_data(new_GuadPhage_samselect)$Treatment, p.adjust.method = "BH")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: erich$Shannon and sample_data(new_GuadPhage_samselect)$Treatment
##
## AA-sucrose AA-blood AA-ZIKV CQ-sucrose CQ-blood
## AA-blood 1.00000 - - - -
## AA-ZIKV 1.00000 1.00000 - - -
## CQ-sucrose 0.00051 0.00051 1.1e-08 - -
## CQ-blood 0.00026 0.00022 5.3e-10 0.08898 -
## CQ-WNV 4.1e-06 4.1e-06 < 2e-16 0.00047 0.14860
##
## P value adjustment method: BH5.3 Phageome of Aedes aegypti
#### color ####
col = plasma(9,alpha = 1, begin=0.9, end = 0, direction = 1)
#### Phyloseq Project ####
GuadPhage_Ae_aegypti <- subset_samples(GuadPhage, mosq_species == "Ae_aegypti" )
GuadPhage_Ae_aegypti## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 263 taxa and 94 samples ]
## sample_data() Sample Data: [ 94 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 263 taxa by 17 taxonomic ranks ]GP.M1_AA = GuadPhage_Ae_aegypti
GP.M1_AA## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 263 taxa and 94 samples ]
## sample_data() Sample Data: [ 94 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 263 taxa by 17 taxonomic ranks ]sum(sample_sums(GP.M1_AA))## [1] 182924sample_data(GP.M1_AA)$Treatment <- factor(sample_data(GP.M1_AA)$Treatment, levels = c("AA-ZIKV","AA-blood","AA-sucrose"))
#### subset ####
sample_sums(GP.M1_AA)## GIM21 GIM22 GIM23 GIM24 GIM25 GIM26 GIM27 GIM28 GIM29 GIM30 GIM31
## 1212 1197 1284 359 1381 1432 1242 1345 1070 1418 1386
## GIM32 GIM33 GIM34 GIM35 GIM36 GIM37 GIM38 GIM39 GIM40 GIM81 GIM82
## 1633 1047 602 1881 1554 1280 1405 1067 1382 1899 996
## GIM83 GIM84 GIM85 GIM86 GIM87 GIM88 GIM89 GIM90 GIM91 GIM92 GIM93
## 330 467 1048 719 1043 438 175 1841 927 0 981
## GIM94 GIM95 GIM96 GIM97 GIM98 GIM99 GIM100 GIM572 GIM573 GIM574 GIM575
## 897 823 1723 1425 1426 817 1418 1317 2604 1496 971
## GIM576 GIM577 GIM578 GIM580 GIM581 GIM582 GIM583 GIM585 GIM586 GIM608 GIM609
## 2124 1525 2095 854 1038 0 2319 1261 1389 1470 598
## GIM610 GIM611 GIM612 GIM613 GIM614 GIM615 GIM616 GIM617 GIM618 GIM619 GIM620
## 1299 1229 434 1062 1176 1251 525 0 2383 1462 586
## GIM621 GIM622 GIM623 GIM624 GIM625 GIM626 GIM627 GIM628 GIM317 GIM318 GIM319
## 848 919 909 666 841 1550 1432 410 732 539 1637
## GIM320 GIM347 GIM349 GIM350 GIM549 GIM553 GIM554 GIM257 GIM258 GIM259 GIM260
## 1688 1339 1088 77045 2418 1118 1259 654 2376 722 456
## GIM284 GIM287 GIM288 GIM290 GIM543 GIM546
## 1188 1686 483 923 150 810sample_sums(GP.M1_AA)[which(sample_sums(GP.M1_AA) == 0)]## GIM92 GIM582 GIM617
## 0 0 0GP.M1_AA_samselect = prune_samples(sample_sums(GP.M1_AA) > 0 , GP.M1_AA)
GP.M1_AA_samselect## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 263 taxa and 91 samples ]
## sample_data() Sample Data: [ 91 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 263 taxa by 17 taxonomic ranks ]GP.M1_AA_select_0 = prune_taxa(taxa_sums(GP.M1_AA_samselect) > 0, GP.M1_AA_samselect)
GP.M1_AA_select_0## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 151 taxa and 91 samples ]
## sample_data() Sample Data: [ 91 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 151 taxa by 17 taxonomic ranks ]#### PCoA ####
GP.ord.GP.M1_AA_samselect <- ordinate(GP.M1_AA_samselect, "PCoA", "bray")
plot_ordination(GP.M1_AA_samselect, GP.ord.GP.M1_AA_samselect, type = "samples", col = "Treatment")+#,label= "name") +
theme_bw() +
geom_point(size = 4, alpha=0.7) +
scale_color_manual(values = c(col[9],col[5],col[2]))+
theme(legend.title=element_blank(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of phages in all Ae. aegypti samples")
set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_AA_samselect, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_AA_samselect))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)##
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 2 0.7636 0.38180 2.2902 0.04947 0.027 *
## Residuals 88 14.6707 0.16671 0.95053
## Total 90 15.4343 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.4 Phageome Culex quinquefasciatus
#### Phyloseq Project ####
GuadPhage_Cx_quinquefasciatus <- subset_samples(GuadPhage, mosq_species == "Cx_quinquefasciatus" )
GuadPhage_Cx_quinquefasciatus## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 263 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 263 taxa by 17 taxonomic ranks ]GP.M1_CQ = GuadPhage_Cx_quinquefasciatus
GP.M1_CQ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 263 taxa and 60 samples ]
## sample_data() Sample Data: [ 60 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 263 taxa by 17 taxonomic ranks ]sample_data(GP.M1_CQ)$group.plaque.qPCR.0.body.0. <- factor(sample_data(GP.M1_CQ)$group.plaque.qPCR.0.body.0.,
levels = c("7dpi-WNV(+/)", "7dpi-WNV(-/+)", "7dpi-WNV(-/-)","CQ-7dpi-blood","CQ-7dpi-sucrose",
"14dpi-WNV(+/)","14dpi-WNV(-/+)","14dpi-WNV(-/-)","CQ-14dpi-blood","CQ-14dpi-sucrose"))
sample_data(GP.M1_CQ)$Info1 <- factor(sample_data(GP.M1_CQ)$Info1, levels = c("CQ-7dpi-WNV", "CQ-7dpi-blood", "CQ-7dpi-sucrose", "CQ-14dpi-WNV", "CQ-14dpi-blood", "CQ-14dpi-sucrose"))
sample_data(GP.M1_CQ)$Treatment <- factor(sample_data(GP.M1_CQ)$Treatment, levels = c("CQ-WNV","CQ-blood","CQ-sucrose"))
levels(sample_data(GP.M1_CQ)$Info1)## [1] "CQ-7dpi-WNV" "CQ-7dpi-blood" "CQ-7dpi-sucrose" "CQ-14dpi-WNV"
## [5] "CQ-14dpi-blood" "CQ-14dpi-sucrose"levels(sample_data(GP.M1_CQ)$group.plaque.qPCR.0.body.0.)## [1] "7dpi-WNV(+/)" "7dpi-WNV(-/+)" "7dpi-WNV(-/-)" "CQ-7dpi-blood"
## [5] "CQ-7dpi-sucrose" "14dpi-WNV(+/)" "14dpi-WNV(-/+)" "14dpi-WNV(-/-)"
## [9] "CQ-14dpi-blood" "CQ-14dpi-sucrose"#### subset ####
sample_sums(GP.M1_CQ)## GIM381 GIM382 GIM383 GIM384 GIM385 GIM386 GIM387 GIM388 GIM389 GIM390 GIM391
## 32527 48435 106963 180814 115000 223177 223468 197054 108111 420751 66868
## GIM392 GIM393 GIM394 GIM395 GIM396 GIM397 GIM398 GIM399 GIM400 GIM421 GIM422
## 14533 717060 216103 262323 389183 129958 153894 423605 35159 406227 80616
## GIM423 GIM424 GIM425 GIM426 GIM427 GIM428 GIM429 GIM430 GIM431 GIM432 GIM433
## 363912 114523 135616 83275 71798 177505 59004 83063 97328 84743 130900
## GIM434 GIM435 GIM436 GIM437 GIM438 GIM439 GIM440 GIM491 GIM492 GIM493 GIM494
## 43898 138201 36543 91653 98097 736556 17718 153484 35892 2698 6944
## GIM495 GIM525 GIM527 GIM528 GIM529 GIM530 GIM451 GIM452 GIM453 GIM454 GIM455
## 4822 26017 10736 79978 21879 0 574667 194018 87396 536174 71076
## GIM471 GIM472 GIM473 GIM475 GIM476
## 479370 572072 7832 369564 599136sample_sums(GP.M1_CQ)[which(sample_sums(GP.M1_CQ) == 0)]## GIM530
## 0GP.M1_CQ_samselect = prune_samples(sample_sums(GP.M1_CQ) > 0 , GP.M1_CQ)
GP.M1_CQ_samselect## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 263 taxa and 59 samples ]
## sample_data() Sample Data: [ 59 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 263 taxa by 17 taxonomic ranks ]GP.M1_CQ_select_0 = prune_taxa(taxa_sums(GP.M1_CQ_samselect) > 0, GP.M1_CQ_samselect)
GP.M1_CQ_select_0## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 164 taxa and 59 samples ]
## sample_data() Sample Data: [ 59 samples by 18 sample variables ]
## tax_table() Taxonomy Table: [ 164 taxa by 17 taxonomic ranks ]#### PCoA all samples ####
GP.ord.GP.M1_CQ_samselect <- ordinate(GP.M1_CQ_samselect, "PCoA", "bray")
plot_ordination(GP.M1_CQ_samselect, GP.ord.GP.M1_CQ_samselect, type = "samples", col = "Treatment")+
theme_bw() +
geom_point(size = 4, alpha=0.7) +
scale_color_manual(values = c(col[9],col[5],col[2]))+
theme(legend.title=element_blank(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of phages in all Cx. quinquefasciatus samples")+
stat_ellipse(type = "norm", linetype = 2,level=0.6,aes(group = Treatment))
set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_CQ_select_0, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_CQ_select_0))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)##
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 2 2.5667 1.28336 6.0101 0.17671 0.001 ***
## Residuals 56 11.9580 0.21353 0.82329
## Total 58 14.5247 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# # Call:
beta <- betadisper(GP.M1_species_select_bray, GP.M1_species_select_sampledf$Treatment)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.0167 0.008348 0.2089 999 0.783
## Residuals 56 2.2379 0.039962#### PCoA 7 dpe ####
GP.M1_CQ_select_0_7dpi = subset_samples(GP.M1_CQ_samselect, dpi == "7dpi")
GP.M1_CQ_select_7dpi = prune_taxa(taxa_sums(GP.M1_CQ_select_0_7dpi) > 0, GP.M1_CQ_select_0_7dpi)
levels(sample_data(GP.M1_CQ_select_7dpi)$group.plaque.qPCR.0.body.0.)## [1] "7dpi-WNV(+/)" "7dpi-WNV(-/+)" "7dpi-WNV(-/-)" "CQ-7dpi-blood"
## [5] "CQ-7dpi-sucrose"GP.ord.GP.M1_CQ_select_7dpi <- ordinate(GP.M1_CQ_select_7dpi, "PCoA", "bray")
plot_ordination(GP.M1_CQ_select_7dpi, GP.ord.GP.M1_CQ_select_7dpi, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.0.body.0.") +
theme_bw()+
geom_point(size = 4, alpha=0.7) +
scale_shape_manual(values=c(shap[1:4],shap[6])) +
scale_color_manual(values = c(col[9],col[5],col[2]))+
theme(legend.title=element_blank(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of phages in Cx. quinquefasciatus at 7 dpe")+
stat_ellipse(type = "norm", linetype = 2,level=0.6,aes(group = Treatment), color="black")
set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_CQ_select_0_7dpi, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_CQ_select_0_7dpi))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)##
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 2 1.2973 0.64865 2.8944 0.17655 0.005 **
## Residuals 27 6.0509 0.22411 0.82345
## Total 29 7.3482 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### PCoA 14 dpe ####
GP.M1_CQ_select_0_14dpi = subset_samples(GP.M1_CQ_samselect, dpi == "14dpi")
GP.M1_CQ_select_14dpi = prune_taxa(taxa_sums(GP.M1_CQ_select_0_14dpi) > 0, GP.M1_CQ_select_0_14dpi)
levels(sample_data(GP.M1_CQ_select_14dpi)$group.plaque.qPCR.0.body.0.)## [1] "14dpi-WNV(+/)" "14dpi-WNV(-/+)" "14dpi-WNV(-/-)" "CQ-14dpi-blood"
## [5] "CQ-14dpi-sucrose"GP.ord.GP.M1_CQ_select_14dpi <- ordinate(GP.M1_CQ_select_14dpi, "PCoA", "bray")
plot_ordination(GP.M1_CQ_select_14dpi, GP.ord.GP.M1_CQ_select_14dpi, type = "samples", color = "Treatment", shape = "group.plaque.qPCR.0.body.0.")+
theme_bw()+
geom_point(size = 4, alpha=0.7) +
scale_shape_manual(values=c(shap[1:4],shap[6])) +
scale_color_manual(values = c(col[9],col[5],col[2]))+
theme(legend.title=element_blank(),
legend.position = "right",
plot.title = element_text(size = 12,hjust = 0.5))+
ggtitle("PCoA of phages in Cx. quinquefasciatus at 14dpi")+
stat_ellipse(type = "norm", linetype = 2,level=0.8,aes(group = Treatment), color="black")
set.seed(1)
# Calculate bray curtis distance matrix
GP.M1_species_select_bray <- phyloseq::distance(GP.M1_CQ_select_0_14dpi, method = "bray")
# make a data frame from the sample_data
GP.M1_species_select_sampledf <- data.frame(sample_data(GP.M1_CQ_select_0_14dpi))
# Adonis test
adonis(GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)##
## Call:
## adonis(formula = GP.M1_species_select_bray ~ Treatment, data = GP.M1_species_select_sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 2 1.9906 0.99531 5.2418 0.28735 0.001 ***
## Residuals 26 4.9368 0.18988 0.71265
## Total 28 6.9275 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.5 Phage correlate with bacteria
#### phage ####
Phage_host = readRDS("~/Desktop/mosq_Guad_infection/RDS/GIM_phage_host.RDS")
Phage_host_subset <- subset_samples(Phage_host, host_genus != "no result")
pOTU <- as.data.frame(t(Phage_host_subset@otu_table))
pTAX <- data.frame(sample_data(Phage_host_subset))
pTAX <- pTAX[,15:16]
pTAX$phageID <- rownames(pTAX)
pOTU$phageID <- rownames(pOTU)
phage_merge <- merge(pOTU, pTAX, by = "phageID")
#### bacteria ####
Bac = readRDS("~/Desktop/mosq_Guad_infection/16s_rRNA/GIM_Bacteriome_decont_SLV132.rds")
BacGenus <- Bac %>%
tax_glom(taxrank = "Genus")
BacOTU <- as.data.frame(t(BacGenus@otu_table))
colnames(BacOTU) <- paste("Bac_", colnames(BacOTU), sep = "")
BacTAX <- as.data.frame(BacGenus@tax_table)
BacOTU$bacID <- rownames(BacOTU)
BacTAX$bacID <- rownames(BacTAX)
Bac_merge <- merge(BacOTU, BacTAX[,5:7], by = "bacID")
Bac_merge$Genus <- gsub("g_", "", Bac_merge$Genus)
#### merge phage + bacteria ####
colnames(phage_merge) <- gsub("host_genus", "Genus", colnames(phage_merge))
Bac_merge_slec <- Bac_merge[Bac_merge$Genus %in% c(phage_merge$Genus, "uc_f_Enterobacteriaceae"),]
phage_merge$Genus_og <- phage_merge$Genus
phage_merge$Genus <- gsub("Enterobacter", "uc_f_Enterobacteriaceae", phage_merge$Genus)
phage_merge$Genus <- gsub("Phytobacter", "uc_f_Enterobacteriaceae", phage_merge$Genus)
phage_Bac <- merge(phage_merge, Bac_merge_slec, all = T, by ="Genus")
phage_Bac <- phage_Bac[, !colnames(phage_Bac) %in% c("Family", "host_species", "bacID", "Bac_NC-p07", "Bac_NC-p06", "Bac_GIM-NC1", "Bac_GIM-NC2")]
phage_Bac_melt <- melt(phage_Bac)## Using Genus, phageID, Genus_og as id variablesphage_Bac_melt <- phage_Bac_melt %>% separate(variable, c("cat","Sample"), sep ="GIM"); head(phage_Bac_melt)## Genus phageID Genus_og cat Sample value
## 1 Asaia GIM439_NODE_511_length_3345_cov_12.118421 Asaia 21 0
## 2 Asaia GIM439_NODE_40_length_9507_cov_18.950795 Asaia 21 0
## 3 Serratia GIM455_NODE_1_length_46241_cov_41.282255 Serratia 21 0
## 4 Serratia GIM476_NODE_105_length_6055_cov_20.624791 Serratia 21 0
## 5 Serratia GIM476_NODE_133_length_5420_cov_30.978664 Serratia 21 0
## 6 Serratia GIM476_NODE_406_length_3177_cov_20.609677 Serratia 21 0phage_Bac_melt$Sample <- gsub("-", "", phage_Bac_melt$Sample)
phage_Bac_melt$Sample <- paste("GIM", phage_Bac_melt$Sample, sep = "")
phage_Bac_melt$cat <- ifelse(phage_Bac_melt$cat == "", "phage", "bacteria"); head(phage_Bac_melt)## Genus phageID Genus_og cat Sample
## 1 Asaia GIM439_NODE_511_length_3345_cov_12.118421 Asaia phage GIM21
## 2 Asaia GIM439_NODE_40_length_9507_cov_18.950795 Asaia phage GIM21
## 3 Serratia GIM455_NODE_1_length_46241_cov_41.282255 Serratia phage GIM21
## 4 Serratia GIM476_NODE_105_length_6055_cov_20.624791 Serratia phage GIM21
## 5 Serratia GIM476_NODE_133_length_5420_cov_30.978664 Serratia phage GIM21
## 6 Serratia GIM476_NODE_406_length_3177_cov_20.609677 Serratia phage GIM21
## value
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0phage_Bac_melt_wide <- dcast(phage_Bac_melt, phageID + Sample + Genus + Genus_og ~ cat, value.var = 'value')
#### add meta ####
head(phage_Bac_melt)## Genus phageID Genus_og cat Sample
## 1 Asaia GIM439_NODE_511_length_3345_cov_12.118421 Asaia phage GIM21
## 2 Asaia GIM439_NODE_40_length_9507_cov_18.950795 Asaia phage GIM21
## 3 Serratia GIM455_NODE_1_length_46241_cov_41.282255 Serratia phage GIM21
## 4 Serratia GIM476_NODE_105_length_6055_cov_20.624791 Serratia phage GIM21
## 5 Serratia GIM476_NODE_133_length_5420_cov_30.978664 Serratia phage GIM21
## 6 Serratia GIM476_NODE_406_length_3177_cov_20.609677 Serratia phage GIM21
## value
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0phage_Bac_melt_wide$Genus <- as.factor(phage_Bac_melt_wide$Genus)
phage_Bac_melt_wide$log10bacteria <- log10(phage_Bac_melt_wide$bacteria)
is.na(phage_Bac_melt_wide$log10bacteria) <- sapply(phage_Bac_melt_wide$log10bacteria, is.infinite)
phage_Bac_melt_wide$log10bacteria[is.na(phage_Bac_melt_wide$log10bacteria)]<-0
phage_Bac_melt_wide$log10phage <- log10(phage_Bac_melt_wide$phage)
is.na(phage_Bac_melt_wide$log10phage) <- sapply(phage_Bac_melt_wide$log10phage, is.infinite)
phage_Bac_melt_wide$log10phage[is.na(phage_Bac_melt_wide$log10phage)]<-0
samplemeta <- data.frame(BacGenus@sam_data); colnames(samplemeta)## [1] "mosq" "body"
## [3] "mosq_or" "head"
## [5] "Index" "mosq_species"
## [7] "Treatment" "dpi"
## [9] "group.plaque.qPCR.1000." "group.plaque.qPCR.0."
## [11] "group.plaque.qPCR.1000.body.1000." "group.plaque.qPCR.0.body.0."
## [13] "Info1" "plaque_assay"
## [15] "qPCR_head._ZIKV_WNV" "qPCR_body._ZIKV_WNV"
## [17] "head_ZIKV.1000" "body_ZIKV.1000"
## [19] "Sample_or_Control" "conc"samplemeta <- samplemeta[ ,colnames(samplemeta) %in% c("mosq", "Index", "mosq_species", "Treatment", "dpi", "plaque_assay", "Info1")]
samplemeta$Sample <- gsub("-", "", rownames(samplemeta))
phage_Bac_melt_wide_meta <- merge(phage_Bac_melt_wide, samplemeta, by= "Sample")
levels(phage_Bac_melt_wide_meta$Genus)## [1] "Asaia" "Serratia"
## [3] "uc_f_Enterobacteriaceae" "Wolbachia"phage_Bac_melt_wide_metaSle <- phage_Bac_melt_wide_meta[phage_Bac_melt_wide_meta$Genus %in% c("Wolbachia", "Serratia"),]
phage_Bac_melt_wide_metaSle$Genus <- factor(phage_Bac_melt_wide_metaSle$Genus, levels = c("Wolbachia", "Serratia"))
phage_Bac_melt_wide_metaSle$mosq_species <- factor(phage_Bac_melt_wide_metaSle$mosq_species, levels = c("Ae_aegypti", "Cx_quinquefasciatus"))#### plot ####
ggscatter(phage_Bac_melt_wide_metaSle, x = "log10bacteria", y = "log10phage",
color = "mosq_species", palette = c("#7294D4", "#C27D38"), size = 2.2, alpha = 0.8,
conf.int = T, cor.coef = TRUE, cor.method = "pearson")+
facet_wrap(~Genus)+
geom_smooth(color="black")+
theme_bw()+
theme(axis.title.y = element_text(size = 12, vjust=0.5),
axis.title.x = element_text(size = 12, hjust=0.5),
panel.grid.minor = element_blank(),
legend.position = "top",
strip.text.x =element_text(size = 14, hjust = 0.5)) +
xlab("log10 host bacteria reads")+
ylab("1og10 phage reads")## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'